From b3d4c8d91722e73e6d66ab05f2ce4e1df828a70c Mon Sep 17 00:00:00 2001 From: misi9170 Date: Mon, 20 May 2024 10:16:58 -0600 Subject: [PATCH 01/37] Adding preliminary tests and some infrastructure for bulk wd heterogeneity. --- floris/core/flow_field.py | 23 ++++++++-- tests/bulk_wd_hetergoeneity_unit_test.py | 54 ++++++++++++++++++++++++ 2 files changed, 74 insertions(+), 3 deletions(-) create mode 100644 tests/bulk_wd_hetergoeneity_unit_test.py diff --git a/floris/core/flow_field.py b/floris/core/flow_field.py index 19b488518..d7ea19616 100644 --- a/floris/core/flow_field.py +++ b/floris/core/flow_field.py @@ -105,8 +105,10 @@ def heterogeneous_config_validator(self, instance: attrs.Attribute, value: dict return # Check that the correct keys are supplied for the heterogeneous_inflow_config dict - for k in ["speed_multipliers", "x", "y"]: - if k not in value.keys(): + speed_het_keys = ["speed_multipliers", "x", "y", "z"] + bulk_wd_het_keys = ["bulk_wd_change", "bulk_wd_x"] + for k in value.keys(): + if k not in speed_het_keys and k not in bulk_wd_het_keys: raise ValueError( "heterogeneous_inflow_config must contain entries for 'speed_multipliers'," f"'x', and 'y', with 'z' optional. Missing '{k}'." @@ -130,7 +132,8 @@ def het_map_validator(self, instance: attrs.Attribute, value: list | None) -> No def __attrs_post_init__(self) -> None: - if self.heterogeneous_inflow_config is not None: + if self.heterogeneous_inflow_config is not None and self._speed_heterogeneity: + print(self._speed_heterogeneity) self.generate_heterogeneous_wind_map() @@ -305,6 +308,20 @@ def generate_heterogeneous_wind_map(self): self.het_map = in_region + @property + def _speed_heterogeneity(self): + if self.heterogeneous_inflow_config is not None: + return "speed_multipliers" in self.heterogeneous_inflow_config + else: + return False + + @property + def _bulk_wd_heterogeneity(self): + if self.heterogeneous_inflow_config is not None: + return "bulk_wd_change" in self.heterogeneous_inflow_config + else: + return False + @staticmethod def interpolate_multiplier_xy(x: NDArrayFloat, y: NDArrayFloat, diff --git a/tests/bulk_wd_hetergoeneity_unit_test.py b/tests/bulk_wd_hetergoeneity_unit_test.py new file mode 100644 index 000000000..ca0b265a5 --- /dev/null +++ b/tests/bulk_wd_hetergoeneity_unit_test.py @@ -0,0 +1,54 @@ +import logging +from pathlib import Path + +import numpy as np +import pytest + +from floris import ( + FlorisModel, + TimeSeries, + WindRose, +) +from floris.optimization.layout_optimization.layout_optimization_base import ( + LayoutOptimization, +) +from floris.optimization.layout_optimization.layout_optimization_scipy import ( + LayoutOptimizationScipy, +) +from floris.wind_data import WindDataBase + + +TEST_DATA = Path(__file__).resolve().parent / "data" +YAML_INPUT = TEST_DATA / "input_full.yaml" + +def test_base_class(caplog): + # Get a test fi (single turbine at 0,0) + fmodel = FlorisModel(configuration=YAML_INPUT) + + # Directly downstream at 270 degrees + sample_x = [500.0] + sample_y = [0.0] + sample_z = [90.0] + + # Sweep across wind directions + wd_array = np.arange(180, 360, 1) + ws_array = 8.0 * np.ones_like(wd_array) + ti_array = 0.06 * np.ones_like(wd_array) + fmodel.set(wind_directions=wd_array, wind_speeds=ws_array, turbulence_intensities=ti_array) + + # Standard case; expect minimum to be at 270 degrees + u_at_points = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z) + assert (wd_array[np.argmin(u_at_points)] == 270) + + # Now, apply bulk wind direction heterogeneity to the flow + fmodel.set( + heterogeneous_inflow_config={ + "bulk_wd_change": [0.0, 10.0], # 10 degree change, CW + "bulk_wd_x": [0, 500.0] + } + ) + + u_at_points = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z) + assert (wd_array[np.argmin(u_at_points)] != 270) + + From 20edf05f79d51ceb5d38c86d0d1d6d6c816d7dd9 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Mon, 20 May 2024 20:47:53 -0600 Subject: [PATCH 02/37] WIP --- floris/core/core.py | 5 +++++ floris/core/grid.py | 25 +++++++++++++++++++++++ floris/core/solver.py | 1 + floris/utilities.py | 26 +++++++++++++++++++++++- tests/bulk_wd_hetergoeneity_unit_test.py | 9 ++++---- 5 files changed, 61 insertions(+), 5 deletions(-) diff --git a/floris/core/core.py b/floris/core/core.py index 89af93bcf..76276bd12 100644 --- a/floris/core/core.py +++ b/floris/core/core.py @@ -96,6 +96,7 @@ def __attrs_post_init__(self) -> None: turbine_diameters=self.farm.rotor_diameters, wind_directions=self.flow_field.wind_directions, grid_resolution=self.solver["turbine_grid_points"], + heterogeneous_inflow_config=self.flow_field.heterogeneous_inflow_config, ) elif self.solver["type"] == "turbine_cubature_grid": self.grid = TurbineCubatureGrid( @@ -103,6 +104,7 @@ def __attrs_post_init__(self) -> None: turbine_diameters=self.farm.rotor_diameters, wind_directions=self.flow_field.wind_directions, grid_resolution=self.solver["turbine_grid_points"], + heterogeneous_inflow_config=self.flow_field.heterogeneous_inflow_config, ) elif self.solver["type"] == "flow_field_grid": self.grid = FlowFieldGrid( @@ -110,6 +112,7 @@ def __attrs_post_init__(self) -> None: turbine_diameters=self.farm.rotor_diameters, wind_directions=self.flow_field.wind_directions, grid_resolution=self.solver["flow_field_grid_points"], + heterogeneous_inflow_config=self.flow_field.heterogeneous_inflow_config, ) elif self.solver["type"] == "flow_field_planar_grid": self.grid = FlowFieldPlanarGrid( @@ -119,6 +122,7 @@ def __attrs_post_init__(self) -> None: normal_vector=self.solver["normal_vector"], planar_coordinate=self.solver["planar_coordinate"], grid_resolution=self.solver["flow_field_grid_points"], + heterogeneous_inflow_config=self.flow_field.heterogeneous_inflow_config, x1_bounds=self.solver["flow_field_bounds"][0], x2_bounds=self.solver["flow_field_bounds"][1], ) @@ -239,6 +243,7 @@ def solve_for_points(self, x, y, z): turbine_coordinates=self.farm.coordinates, turbine_diameters=self.farm.rotor_diameters, wind_directions=self.flow_field.wind_directions, + heterogeneous_inflow_config=self.flow_field.heterogeneous_inflow_config, grid_resolution=1, x_center_of_rotation=self.grid.x_center_of_rotation, y_center_of_rotation=self.grid.y_center_of_rotation diff --git a/floris/core/grid.py b/floris/core/grid.py index b20ca25ce..46a90dee8 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -18,6 +18,7 @@ from floris.utilities import ( reverse_rotate_coordinates_rel_west, rotate_coordinates_rel_west, + warp_grid_for_wind_direction_heterogeneity, ) @@ -52,6 +53,7 @@ class Grid(ABC, BaseClass): turbine_diameters: NDArrayFloat = field(converter=floris_array_converter) wind_directions: NDArrayFloat = field(converter=floris_array_converter) grid_resolution: int | Iterable = field() + heterogeneous_inflow_config: dict | None = field() n_turbines: int = field(init=False) n_findex: int = field(init=False) @@ -187,6 +189,29 @@ def set_grid(self) -> None: self.turbine_coordinates, ) + if (self.heterogeneous_inflow_config is not None + and "bulk_wd_change" in self.heterogeneous_inflow_config): + wd_het_x_points = np.array(self.heterogeneous_inflow_config["bulk_wd_x"]) + wd_het_values = np.array(self.heterogeneous_inflow_config["bulk_wd_change"]) + # TODO: does the below work as expected? Do I need the y values too? + coordinates_to_rotate = np.concatenate( + ( + wd_het_x_points[:,:,None], + np.zeros((wd_het_x_points.shape[0], wd_het_x_points.shape[1], 2)) + ), + axis=2 + )[None,:,:,:] + wd_het_x_points = rotate_coordinates_rel_west( + np.repeat(self.wind_directions[None,:], wd_het_x_points.shape[1], axis=0), + coordinates_to_rotate, + self.x_center_of_rotation, + self.y_center_of_rotation, + )[0][:, :, 0].T + + x, y, z = warp_grid_for_wind_direction_heterogeneity( + x, y, z, wd_het_x_points, wd_het_values + ) + # - **rloc** (*float, optional): A value, from 0 to 1, that determines # the width/height of the grid of points on the rotor as a ratio of # the rotor radius. diff --git a/floris/core/solver.py b/floris/core/solver.py index a7c3d8796..42479f734 100644 --- a/floris/core/solver.py +++ b/floris/core/solver.py @@ -287,6 +287,7 @@ def full_flow_sequential_solver( turbine_coordinates=turbine_grid_farm.coordinates, turbine_diameters=turbine_grid_farm.rotor_diameters, wind_directions=turbine_grid_flow_field.wind_directions, + heterogeneous_inflow_config=turbine_grid_flow_field.heterogeneous_inflow_config, grid_resolution=3, ) turbine_grid_farm.expand_farm_properties( diff --git a/floris/utilities.py b/floris/utilities.py index 074d9a1b3..89e53289d 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -139,7 +139,7 @@ def rotate_coordinates_rel_west( # Calculate the difference in given wind direction from 270 / West wind_deviation_from_west = wind_delta(wind_directions) - wind_deviation_from_west = np.reshape(wind_deviation_from_west, (len(wind_directions), 1)) + wind_deviation_from_west = np.expand_dims(wind_deviation_from_west, axis=-1) # Construct the arrays storing the turbine locations x_coordinates, y_coordinates, z_coordinates = coordinates.T @@ -224,6 +224,30 @@ def reverse_rotate_coordinates_rel_west( return grid_x_reversed, grid_y_reversed, grid_z_reversed +def warp_grid_for_wind_direction_heterogeneity( + x, + y, + z, + wd_het_x_points, + wd_het_values, + ): + + x0 = wd_het_x_points[:,0:1] + + # Compute radius of curvature for the warp (not sign swap in tand) + import ipdb; ipdb.set_trace() + warp_radius = (wd_het_x_points[:,1] - x0) / tand(wd_het_values[:,0] - wd_het_values[:,1]) + + # Compute the angle of warping for each grid point (radians) + warp_theta = np.arctan2(x - x0, warp_radius) + + # Apply warp + x_warped = warp_radius * np.sin(warp_theta) + y * np.sin(warp_theta) + x0 + y_warped = warp_radius * np.cos(warp_theta) + y * np.cos(warp_theta) - warp_radius + + # Return warped grid + return np.where(x >= x0, x_warped, x), np.where(x >= x0, y_warped, y), z + class Loader(yaml.SafeLoader): diff --git a/tests/bulk_wd_hetergoeneity_unit_test.py b/tests/bulk_wd_hetergoeneity_unit_test.py index ca0b265a5..9f69d4d7e 100644 --- a/tests/bulk_wd_hetergoeneity_unit_test.py +++ b/tests/bulk_wd_hetergoeneity_unit_test.py @@ -24,6 +24,7 @@ def test_base_class(caplog): # Get a test fi (single turbine at 0,0) fmodel = FlorisModel(configuration=YAML_INPUT) + fmodel.set(layout_x=[0, 1000], layout_y=[0, 0]) # Directly downstream at 270 degrees sample_x = [500.0] @@ -43,12 +44,12 @@ def test_base_class(caplog): # Now, apply bulk wind direction heterogeneity to the flow fmodel.set( heterogeneous_inflow_config={ - "bulk_wd_change": [0.0, 10.0], # 10 degree change, CW - "bulk_wd_x": [0, 500.0] + "bulk_wd_change": [[0.0, -10.0]]*180, # -10 degree change, CW + "bulk_wd_x": [[0, 500.0]]*180 } - ) + ) # TODO: Build something that checks the dimensions of the inputs here u_at_points = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z) - assert (wd_array[np.argmin(u_at_points)] != 270) + assert (wd_array[np.argmin(u_at_points)] == 280) From d73eb29599fdd4deb462245206eee13b0ebf6600 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Tue, 21 May 2024 10:18:40 -0600 Subject: [PATCH 03/37] Adding test for turbines only. --- floris/core/grid.py | 2 +- floris/utilities.py | 3 +-- tests/bulk_wd_hetergoeneity_unit_test.py | 31 ++++++++++++++++++++++-- 3 files changed, 31 insertions(+), 5 deletions(-) diff --git a/floris/core/grid.py b/floris/core/grid.py index 46a90dee8..b7230d10c 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -189,7 +189,7 @@ def set_grid(self) -> None: self.turbine_coordinates, ) - if (self.heterogeneous_inflow_config is not None + if (self.heterogeneous_inflow_config is not None and "bulk_wd_change" in self.heterogeneous_inflow_config): wd_het_x_points = np.array(self.heterogeneous_inflow_config["bulk_wd_x"]) wd_het_values = np.array(self.heterogeneous_inflow_config["bulk_wd_change"]) diff --git a/floris/utilities.py b/floris/utilities.py index 89e53289d..406d1b91c 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -235,8 +235,7 @@ def warp_grid_for_wind_direction_heterogeneity( x0 = wd_het_x_points[:,0:1] # Compute radius of curvature for the warp (not sign swap in tand) - import ipdb; ipdb.set_trace() - warp_radius = (wd_het_x_points[:,1] - x0) / tand(wd_het_values[:,0] - wd_het_values[:,1]) + warp_radius = (wd_het_x_points[:,1:] - x0) / tand(wd_het_values[:,0:1] - wd_het_values[:,1:]) # Compute the angle of warping for each grid point (radians) warp_theta = np.arctan2(x - x0, warp_radius) diff --git a/tests/bulk_wd_hetergoeneity_unit_test.py b/tests/bulk_wd_hetergoeneity_unit_test.py index 9f69d4d7e..80134ead8 100644 --- a/tests/bulk_wd_hetergoeneity_unit_test.py +++ b/tests/bulk_wd_hetergoeneity_unit_test.py @@ -21,10 +21,9 @@ TEST_DATA = Path(__file__).resolve().parent / "data" YAML_INPUT = TEST_DATA / "input_full.yaml" -def test_base_class(caplog): +def test_bulk_wd_heterogeneity_flow_field(): # Get a test fi (single turbine at 0,0) fmodel = FlorisModel(configuration=YAML_INPUT) - fmodel.set(layout_x=[0, 1000], layout_y=[0, 0]) # Directly downstream at 270 degrees sample_x = [500.0] @@ -52,4 +51,32 @@ def test_base_class(caplog): u_at_points = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z) assert (wd_array[np.argmin(u_at_points)] == 280) +def test_bulk_wd_heterogeneity_turbine(): + fmodel = FlorisModel(configuration=YAML_INPUT) + fmodel.set(layout_x=[0, 600], layout_y=[0, 0]) + + # Run straight as a benchmark + fmodel.set(wind_directions=[270.0], wind_speeds=[8.0], turbulence_intensities=[0.06]) + fmodel.run() + powers_straight = fmodel.get_turbine_powers()[0,:] + + + # Single wind direction, with two wind direction shifts as well as straight + fmodel.set( + wind_directions=[270.0, 270.0, 270.0], + wind_speeds=[8.0, 8.0, 8.0], + turbulence_intensities=[0.06, 0.06, 0.06], + heterogeneous_inflow_config={ + "bulk_wd_change": [[0.0, 0.0], [0.0, -10.0], [0.0, 10.0]], # 10 degree changes + "bulk_wd_x": [[0, 1000], [0, 1000], [0, 1000]] # Past downstream turbine + } + ) + + fmodel.run() + powers = fmodel.get_turbine_powers() + + assert (powers_straight == powers[0,:]).all() # Verify straight case + assert (powers[:,0] == powers[0,0]).all() # Upstream turbine not affected + assert powers[0,1] < powers[1,1] # wake shifted away from downstream turbine + assert powers[1,1] == powers[2,1] # Shifted wake is symmetric From 32701358e9de601cff38fe6cf73cacc6264b4d31 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Tue, 21 May 2024 12:33:12 -0600 Subject: [PATCH 04/37] switched to sliding method. --- floris/core/grid.py | 29 +++++----- floris/utilities.py | 36 +++++++++--- tests/bulk_wd_hetergoeneity_unit_test.py | 71 ++++++++++++++---------- 3 files changed, 82 insertions(+), 54 deletions(-) diff --git a/floris/core/grid.py b/floris/core/grid.py index b7230d10c..00577a35e 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -194,22 +194,21 @@ def set_grid(self) -> None: wd_het_x_points = np.array(self.heterogeneous_inflow_config["bulk_wd_x"]) wd_het_values = np.array(self.heterogeneous_inflow_config["bulk_wd_change"]) # TODO: does the below work as expected? Do I need the y values too? - coordinates_to_rotate = np.concatenate( - ( - wd_het_x_points[:,:,None], - np.zeros((wd_het_x_points.shape[0], wd_het_x_points.shape[1], 2)) - ), - axis=2 - )[None,:,:,:] - wd_het_x_points = rotate_coordinates_rel_west( - np.repeat(self.wind_directions[None,:], wd_het_x_points.shape[1], axis=0), - coordinates_to_rotate, - self.x_center_of_rotation, - self.y_center_of_rotation, - )[0][:, :, 0].T - + # coordinates_to_rotate = np.concatenate( + # ( + # wd_het_x_points[:,:,None], + # np.zeros((wd_het_x_points.shape[0], wd_het_x_points.shape[1], 2)) + # ), + # axis=2 + # )[None,:,:,:] + # wd_het_x_points = rotate_coordinates_rel_west( + # np.repeat(self.wind_directions[None,:], wd_het_x_points.shape[1], axis=0), + # coordinates_to_rotate, + # self.x_center_of_rotation, + # self.y_center_of_rotation, + # )[0][:, :, 0].T x, y, z = warp_grid_for_wind_direction_heterogeneity( - x, y, z, wd_het_x_points, wd_het_values + x, y, z, wd_het_x_points + x.min(), wd_het_values ) # - **rloc** (*float, optional): A value, from 0 to 1, that determines diff --git a/floris/utilities.py b/floris/utilities.py index 406d1b91c..1aecfe7cc 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -233,19 +233,37 @@ def warp_grid_for_wind_direction_heterogeneity( ): x0 = wd_het_x_points[:,0:1] + x1 = wd_het_x_points[:,1:] # Compute radius of curvature for the warp (not sign swap in tand) - warp_radius = (wd_het_x_points[:,1:] - x0) / tand(wd_het_values[:,0:1] - wd_het_values[:,1:]) - - # Compute the angle of warping for each grid point (radians) - warp_theta = np.arctan2(x - x0, warp_radius) - - # Apply warp - x_warped = warp_radius * np.sin(warp_theta) + y * np.sin(warp_theta) + x0 - y_warped = warp_radius * np.cos(warp_theta) + y * np.cos(warp_theta) - warp_radius + warp_radius = (x1 - x0) / tand(wd_het_values[:,0:1] - wd_het_values[:,1:]) + + ### Below is code for a "true" warp in x and y, but is not complete. + # # Compute the angle of warping for each grid point (radians) + # warp_theta = np.arctan2(x - x0, np.abs(warp_radius))*np.sign(warp_radius) + + # # Apply warp + # import ipdb; ipdb.set_trace() + # x_warped = np.where( # Also add end points + # (x >= x0) & (x <= x1) & (np.abs(warp_radius) < np.inf), + # warp_radius * np.sin(warp_theta) + y * np.sin(warp_theta) + x0, + # x + # ) + # y_warped = np.where( + # (x >= x0) & (x <= x1) & (np.abs(warp_radius) < np.inf), + # warp_radius * np.cos(warp_theta) + y * np.cos(warp_theta) - warp_radius, + # y + # ) + + ### Instead, use a simpler y "slide" approach + y_warped = np.where( + (np.abs(warp_radius) != np.inf ) & (x >= x0) & (x <= x1), + warp_radius - np.sqrt(warp_radius**2 - (x - x0)**2) + y, + y + ) # Return warped grid - return np.where(x >= x0, x_warped, x), np.where(x >= x0, y_warped, y), z + return x, y_warped, z class Loader(yaml.SafeLoader): diff --git a/tests/bulk_wd_hetergoeneity_unit_test.py b/tests/bulk_wd_hetergoeneity_unit_test.py index 80134ead8..3d8721062 100644 --- a/tests/bulk_wd_hetergoeneity_unit_test.py +++ b/tests/bulk_wd_hetergoeneity_unit_test.py @@ -21,6 +21,47 @@ TEST_DATA = Path(__file__).resolve().parent / "data" YAML_INPUT = TEST_DATA / "input_full.yaml" +def test_bulk_wd_heterogeneity_turbine(): + fmodel = FlorisModel(configuration=YAML_INPUT) + fmodel.set(layout_x=[0, 600], layout_y=[0, 0]) + + # Run straight as a benchmark + fmodel.set(wind_directions=[270.0], wind_speeds=[8.0], turbulence_intensities=[0.06]) + fmodel.run() + powers_straight = fmodel.get_turbine_powers()[0,:] + + + # Single wind direction, with two wind direction shifts as well as straight + fmodel.set( + wind_directions=[270.0, 270.0, 270.0], + wind_speeds=[8.0, 8.0, 8.0], + turbulence_intensities=[0.06, 0.06, 0.06], + heterogeneous_inflow_config={ + "bulk_wd_change": [[0.0, 0.0], [0.0, -10.0], [0.0, 10.0]], # 10 degree changes + "bulk_wd_x": [[0, 1000], [0, 1000], [0, 1000]] # Past downstream turbine + } + ) + + fmodel.run() + powers = fmodel.get_turbine_powers() + + assert (powers_straight == powers[0,:]).all() # Verify straight case + + assert (powers[:,0] == powers[0,0]).all() # Upstream turbine not affected + assert powers[0,1] < powers[1,1] # wake shifted away from downstream turbine + assert np.allclose(powers[1,1], powers[2,1]) # Shifted wake is symmetric + + # Check works for other directions too + fmodel.set( + wind_directions=[225.0, 225.0, 225.0], + layout_x=[0, 600/np.sqrt(2)], + layout_y=[0, 600/np.sqrt(2)], + ) + fmodel.run() + powers_2 = fmodel.get_turbine_powers() + + assert np.allclose(powers_2, powers) + def test_bulk_wd_heterogeneity_flow_field(): # Get a test fi (single turbine at 0,0) fmodel = FlorisModel(configuration=YAML_INPUT) @@ -50,33 +91,3 @@ def test_bulk_wd_heterogeneity_flow_field(): u_at_points = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z) assert (wd_array[np.argmin(u_at_points)] == 280) - -def test_bulk_wd_heterogeneity_turbine(): - fmodel = FlorisModel(configuration=YAML_INPUT) - fmodel.set(layout_x=[0, 600], layout_y=[0, 0]) - - # Run straight as a benchmark - fmodel.set(wind_directions=[270.0], wind_speeds=[8.0], turbulence_intensities=[0.06]) - fmodel.run() - powers_straight = fmodel.get_turbine_powers()[0,:] - - - # Single wind direction, with two wind direction shifts as well as straight - fmodel.set( - wind_directions=[270.0, 270.0, 270.0], - wind_speeds=[8.0, 8.0, 8.0], - turbulence_intensities=[0.06, 0.06, 0.06], - heterogeneous_inflow_config={ - "bulk_wd_change": [[0.0, 0.0], [0.0, -10.0], [0.0, 10.0]], # 10 degree changes - "bulk_wd_x": [[0, 1000], [0, 1000], [0, 1000]] # Past downstream turbine - } - ) - - fmodel.run() - powers = fmodel.get_turbine_powers() - - assert (powers_straight == powers[0,:]).all() # Verify straight case - - assert (powers[:,0] == powers[0,0]).all() # Upstream turbine not affected - assert powers[0,1] < powers[1,1] # wake shifted away from downstream turbine - assert powers[1,1] == powers[2,1] # Shifted wake is symmetric From 3912720cf5395537dc36ac63d80d809f5924c560 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Tue, 21 May 2024 13:57:54 -0600 Subject: [PATCH 05/37] Some running visualization examples now. --- examples/examples_heterogeneous/bulk_het.py | 79 +++++++++++++++++++ examples/examples_heterogeneous/bulk_het_2.py | 48 +++++++++++ floris/core/grid.py | 48 ++++++++--- floris/floris_model.py | 24 +++--- 4 files changed, 177 insertions(+), 22 deletions(-) create mode 100644 examples/examples_heterogeneous/bulk_het.py create mode 100644 examples/examples_heterogeneous/bulk_het_2.py diff --git a/examples/examples_heterogeneous/bulk_het.py b/examples/examples_heterogeneous/bulk_het.py new file mode 100644 index 000000000..f2080c76a --- /dev/null +++ b/examples/examples_heterogeneous/bulk_het.py @@ -0,0 +1,79 @@ +import numpy as np +import matplotlib.pyplot as plt +from floris import FlorisModel + +from pathlib import Path + + +# Get a test fi (single turbine at 0,0) +fmodel = FlorisModel("inputs/gch.yaml") +fmodel.set(layout_x=[0], layout_y=[0]) + +# Directly downstream at 270 degrees +option = False +if not option: + sample_x = [500.0] + sample_y = [0.0] + sample_z = [90.0] + + # Sweep across wind directions + wd_array = np.arange(180, 360, 1) + ws_array = 8.0 * np.ones_like(wd_array) + ti_array = 0.06 * np.ones_like(wd_array) + fmodel.set(wind_directions=wd_array, wind_speeds=ws_array, turbulence_intensities=ti_array) + + # Standard case; expect minimum to be at 270 degrees + u_at_points_0 = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z) + + # Now, apply bulk wind direction heterogeneity to the flow + fmodel.set( + heterogeneous_inflow_config={ + "bulk_wd_change": [[0.0, 10.0]]*180, # -10 degree change, CW + "bulk_wd_x": [[0, 500.0]]*180 + } + ) # TODO: Build something that checks the dimensions of the inputs here + u_at_points_1 = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z) + + fmodel.set( + heterogeneous_inflow_config={ + "bulk_wd_change": [[0.0, -10.0]]*180, # -10 degree change, CW + "bulk_wd_x": [[0, 500.0]]*180 + } + ) # TODO: Build something that checks the dimensions of the inputs here + u_at_points_2 = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z) + + + fig, ax = plt.subplots(1,1) + ax.plot(wd_array, u_at_points_0, label="Standard", color="black") + ax.plot(wd_array, u_at_points_1, label="Shifted 10") + ax.plot(wd_array, u_at_points_2, label="Shifted -10") + +else: + ### Let's try something else + sample_x = [500.0]*100 + sample_y = np.linspace(-500, 500, 100) + sample_z = [90.0]*100 + + fmodel.set( + wind_directions=[270.0], + wind_speeds=[8.0], + turbulence_intensities=[0.06], + heterogeneous_inflow_config={} + ) + u_at_points_00 = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z).flatten() + + fmodel.set( + heterogeneous_inflow_config={ + "bulk_wd_change": [[0.0, -30.0]], # -10 degree change, CW + "bulk_wd_x": [[0, 500.0]] + }, + ) + u_at_points_01 = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z).flatten() + + fig, ax = plt.subplots(1,1) + ax.plot(sample_y, u_at_points_00, label="Standard", color="black") + ax.plot(sample_y, u_at_points_01, label="Shifted 10") + +plt.show() + + diff --git a/examples/examples_heterogeneous/bulk_het_2.py b/examples/examples_heterogeneous/bulk_het_2.py new file mode 100644 index 000000000..9f4274026 --- /dev/null +++ b/examples/examples_heterogeneous/bulk_het_2.py @@ -0,0 +1,48 @@ +import numpy as np +import matplotlib.pyplot as plt + +import floris.layout_visualization as layoutviz +from floris import FlorisModel +from floris.flow_visualization import visualize_cut_plane + + +# Get a test fi (single turbine at 0,0) +fmodel = FlorisModel("inputs/gch.yaml") +fmodel.set(layout_x=[0], layout_y=[0]) + +# Directly downstream at 270 degrees +### Let's try something else +sample_x = [500.0]*100 +sample_y = np.linspace(-500, 500, 100) +sample_z = [90.0]*100 + +fmodel.set( + wind_directions=[270.0], + wind_speeds=[8.0], + turbulence_intensities=[0.06], + heterogeneous_inflow_config={ + "bulk_wd_change": [[0.0, 15.0]], # -10 degree change, CW + "bulk_wd_x": [[0, 2000.0]] + }, +) +horizontal_plane = fmodel.calculate_horizontal_plane( + x_resolution=200, + y_resolution=100, + height=90.0, +) + +# Plot the flow field with rotors +fig, ax = plt.subplots() +visualize_cut_plane( + horizontal_plane, + ax=ax, + label_contours=False, + title="Horizontal Flow with Turbine Rotors and labels", +) + +# Plot the turbine rotors +layoutviz.plot_turbine_rotors(fmodel, ax=ax) + +plt.show() + + diff --git a/floris/core/grid.py b/floris/core/grid.py index 00577a35e..fa01c600e 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -316,6 +316,13 @@ def set_grid(self) -> None: self.wind_directions, self.turbine_coordinates ) + if (self.heterogeneous_inflow_config is not None + and "bulk_wd_change" in self.heterogeneous_inflow_config): + wd_het_x_points = np.array(self.heterogeneous_inflow_config["bulk_wd_x"]) + wd_het_values = np.array(self.heterogeneous_inflow_config["bulk_wd_change"]) + x, y, z = warp_grid_for_wind_direction_heterogeneity( + x, y, z, wd_het_x_points + x.min(), wd_het_values + ) # Coefficients cubature_coefficients = TurbineCubatureGrid.get_cubature_coefficients(self.grid_resolution) @@ -487,6 +494,13 @@ def set_grid(self) -> None: self.wind_directions, self.turbine_coordinates ) + if (self.heterogeneous_inflow_config is not None + and "bulk_wd_change" in self.heterogeneous_inflow_config): + wd_het_x_points = np.array(self.heterogeneous_inflow_config["bulk_wd_x"]) + wd_het_values = np.array(self.heterogeneous_inflow_config["bulk_wd_change"]) + x, y, z = warp_grid_for_wind_direction_heterogeneity( + x, y, z, wd_het_x_points + x.min(), wd_het_values + ) # Construct the arrays storing the grid points eps = 0.01 @@ -582,10 +596,6 @@ def set_grid(self) -> None: indexing="ij" ) - self.x_sorted = x_points[None, :, :, :] - self.y_sorted = y_points[None, :, :, :] - self.z_sorted = z_points[None, :, :, :] - elif self.normal_vector == "x": # Rules of thumb for cross plane if self.x1_bounds is None: self.x1_bounds = (np.min(y) - 2 * max_diameter, np.max(y) + 2 * max_diameter) @@ -600,10 +610,6 @@ def set_grid(self) -> None: indexing="ij" ) - self.x_sorted = x_points[None, :, :, :] - self.y_sorted = y_points[None, :, :, :] - self.z_sorted = z_points[None, :, :, :] - elif self.normal_vector == "y": # Rules of thumb for y plane if self.x1_bounds is None: self.x1_bounds = (np.min(x) - 2 * max_diameter, np.max(x) + 10 * max_diameter) @@ -618,9 +624,10 @@ def set_grid(self) -> None: indexing="ij" ) - self.x_sorted = x_points[None, :, :, :] - self.y_sorted = y_points[None, :, :, :] - self.z_sorted = z_points[None, :, :, :] + self.x_sorted = x_points[None, :, :, :] + self.y_sorted = y_points[None, :, :, :] + self.z_sorted = z_points[None, :, :, :] + # Now calculate grid coordinates in original frame (from 270 deg perspective) self.x_sorted_inertial_frame, self.y_sorted_inertial_frame, self.z_sorted_inertial_frame = \ @@ -633,6 +640,18 @@ def set_grid(self) -> None: y_center_of_rotation=self.y_center_of_rotation, ) + if (self.heterogeneous_inflow_config is not None + and "bulk_wd_change" in self.heterogeneous_inflow_config): + wd_het_x_points = np.array(self.heterogeneous_inflow_config["bulk_wd_x"]) + wd_het_values = np.array(self.heterogeneous_inflow_config["bulk_wd_change"]) + x_points, y_points, z_points = warp_grid_for_wind_direction_heterogeneity( + x_points, y_points, z_points, wd_het_x_points + x.min(), wd_het_values + ) + + self.x_sorted = x_points[None, :, :, :] + self.y_sorted = y_points[None, :, :, :] + self.z_sorted = z_points[None, :, :, :] + @define class PointsGrid(Grid): """ @@ -679,6 +698,13 @@ def set_grid(self) -> None: x_center_of_rotation=self.x_center_of_rotation, y_center_of_rotation=self.y_center_of_rotation ) + if (self.heterogeneous_inflow_config is not None + and "bulk_wd_change" in self.heterogeneous_inflow_config): + wd_het_x_points = np.array(self.heterogeneous_inflow_config["bulk_wd_x"]) + wd_het_values = np.array(self.heterogeneous_inflow_config["bulk_wd_change"]) + x, y, z = warp_grid_for_wind_direction_heterogeneity( + x, y, z, wd_het_x_points, wd_het_values + ) self.x_sorted = x[:,:,None,None] self.y_sorted = y[:,:,None,None] self.z_sorted = z[:,:,None,None] diff --git a/floris/floris_model.py b/floris/floris_model.py index 1d05f3d3e..1b4faf032 100644 --- a/floris/floris_model.py +++ b/floris/floris_model.py @@ -957,16 +957,18 @@ def set_for_viz(self, findex: int, solver_settings: dict) -> None: # If not None, set the heterogeneous inflow configuration if self.core.flow_field.heterogeneous_inflow_config is not None: - heterogeneous_inflow_config = { - 'x': self.core.flow_field.heterogeneous_inflow_config['x'], - 'y': self.core.flow_field.heterogeneous_inflow_config['y'], - 'speed_multipliers': - self.core.flow_field.heterogeneous_inflow_config['speed_multipliers'][findex:findex+1], - } - if 'z' in self.core.flow_field.heterogeneous_inflow_config: - heterogeneous_inflow_config['z'] = ( - self.core.flow_field.heterogeneous_inflow_config['z'] - ) + speed_het_keys = ["speed_multipliers", "x", "y", "z"] + bulk_wd_het_keys = ["bulk_wd_change", "bulk_wd_x"] + heterogeneous_inflow_config = {} + for k in self.core.flow_field.heterogeneous_inflow_config.keys(): + if ( + (k in speed_het_keys + or k in bulk_wd_het_keys) + and self.core.flow_field.heterogeneous_inflow_config[k] is not None + ): + heterogeneous_inflow_config[k] = ( + self.core.flow_field.heterogeneous_inflow_config[k][findex:findex+1] + ) else: heterogeneous_inflow_config = None @@ -978,7 +980,7 @@ def set_for_viz(self, findex: int, solver_settings: dict) -> None: power_setpoints=self.core.farm.power_setpoints[findex:findex+1,:], awc_modes=self.core.farm.awc_modes[findex:findex+1,:], awc_amplitudes=self.core.farm.awc_amplitudes[findex:findex+1,:], - heterogeneous_inflow_config = heterogeneous_inflow_config, + heterogeneous_inflow_config=heterogeneous_inflow_config, solver_settings=solver_settings, ) From 3ddb1de188bab450154790b4fefc752e55dcc9f8 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Tue, 21 May 2024 14:48:55 -0600 Subject: [PATCH 06/37] Fix sign issue. --- examples/examples_heterogeneous/bulk_het_2.py | 6 +++--- floris/utilities.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/examples_heterogeneous/bulk_het_2.py b/examples/examples_heterogeneous/bulk_het_2.py index 9f4274026..27c119650 100644 --- a/examples/examples_heterogeneous/bulk_het_2.py +++ b/examples/examples_heterogeneous/bulk_het_2.py @@ -7,8 +7,8 @@ # Get a test fi (single turbine at 0,0) -fmodel = FlorisModel("inputs/gch.yaml") -fmodel.set(layout_x=[0], layout_y=[0]) +fmodel = FlorisModel("../inputs/gch.yaml") +fmodel.set(layout_x=[0, 600], layout_y=[0, 0]) # Directly downstream at 270 degrees ### Let's try something else @@ -21,7 +21,7 @@ wind_speeds=[8.0], turbulence_intensities=[0.06], heterogeneous_inflow_config={ - "bulk_wd_change": [[0.0, 15.0]], # -10 degree change, CW + "bulk_wd_change": [[0.0, -15.0]], # -10 degree change, CW "bulk_wd_x": [[0, 2000.0]] }, ) diff --git a/floris/utilities.py b/floris/utilities.py index 1aecfe7cc..abb14e5b9 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -258,7 +258,7 @@ def warp_grid_for_wind_direction_heterogeneity( ### Instead, use a simpler y "slide" approach y_warped = np.where( (np.abs(warp_radius) != np.inf ) & (x >= x0) & (x <= x1), - warp_radius - np.sqrt(warp_radius**2 - (x - x0)**2) + y, + -warp_radius + np.sign(warp_radius) * np.sqrt(warp_radius**2 - (x - x0)**2) + y, y ) From 6469c6372f07cb794105eff30a50ecc90abd6537 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Thu, 23 May 2024 09:31:46 -0600 Subject: [PATCH 07/37] upddated method. --- examples/examples_heterogeneous/bulk_het_2.py | 9 +-- floris/utilities.py | 71 +++++++++++-------- 2 files changed, 47 insertions(+), 33 deletions(-) diff --git a/examples/examples_heterogeneous/bulk_het_2.py b/examples/examples_heterogeneous/bulk_het_2.py index 27c119650..645e3b9cd 100644 --- a/examples/examples_heterogeneous/bulk_het_2.py +++ b/examples/examples_heterogeneous/bulk_het_2.py @@ -21,8 +21,9 @@ wind_speeds=[8.0], turbulence_intensities=[0.06], heterogeneous_inflow_config={ - "bulk_wd_change": [[0.0, -15.0]], # -10 degree change, CW - "bulk_wd_x": [[0, 2000.0]] + "bulk_wd_change": [[0.0, -20.0, 20.0, -20.0]], # -10 degree change, CW + #"bulk_wd_change": [[0.0, 0.0, 0.0, 0.0]], # -10 degree change, CCW + "bulk_wd_x": [[0, 600, 1200, 1800]] }, ) horizontal_plane = fmodel.calculate_horizontal_plane( @@ -37,11 +38,11 @@ horizontal_plane, ax=ax, label_contours=False, - title="Horizontal Flow with Turbine Rotors and labels", + #title="Horizontal Flow with Turbine Rotors and labels", ) # Plot the turbine rotors -layoutviz.plot_turbine_rotors(fmodel, ax=ax) +ax.grid() plt.show() diff --git a/floris/utilities.py b/floris/utilities.py index abb14e5b9..6bc88587c 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -232,35 +232,48 @@ def warp_grid_for_wind_direction_heterogeneity( wd_het_values, ): - x0 = wd_het_x_points[:,0:1] - x1 = wd_het_x_points[:,1:] - - # Compute radius of curvature for the warp (not sign swap in tand) - warp_radius = (x1 - x0) / tand(wd_het_values[:,0:1] - wd_het_values[:,1:]) - - ### Below is code for a "true" warp in x and y, but is not complete. - # # Compute the angle of warping for each grid point (radians) - # warp_theta = np.arctan2(x - x0, np.abs(warp_radius))*np.sign(warp_radius) - - # # Apply warp - # import ipdb; ipdb.set_trace() - # x_warped = np.where( # Also add end points - # (x >= x0) & (x <= x1) & (np.abs(warp_radius) < np.inf), - # warp_radius * np.sin(warp_theta) + y * np.sin(warp_theta) + x0, - # x - # ) - # y_warped = np.where( - # (x >= x0) & (x <= x1) & (np.abs(warp_radius) < np.inf), - # warp_radius * np.cos(warp_theta) + y * np.cos(warp_theta) - warp_radius, - # y - # ) - - ### Instead, use a simpler y "slide" approach - y_warped = np.where( - (np.abs(warp_radius) != np.inf ) & (x >= x0) & (x <= x1), - -warp_radius + np.sign(warp_radius) * np.sqrt(warp_radius**2 - (x - x0)**2) + y, - y - ) + y_warped = y.copy() + + for i in range(wd_het_values.shape[1]-1): + + x0 = wd_het_x_points[:,i:i+1] + x1 = wd_het_x_points[:,i+1:i+2] + + # Compute radius of curvature for the warp (not sign swap in tand) + warp_radius = (x1 - x0) / tand(wd_het_values[:,i:i+1] - wd_het_values[:,i+1:i+2]) + + ### Below is code for a "true" warp in x and y, but is not complete. + # # Compute the angle of warping for each grid point (radians) + # warp_theta = np.arctan2(x - x0, np.abs(warp_radius))*np.sign(warp_radius) + + # # Apply warp + # import ipdb; ipdb.set_trace() + # x_warped = np.where( # Also add end points + # (x >= x0) & (x <= x1) & (np.abs(warp_radius) < np.inf), + # warp_radius * np.sin(warp_theta) + y * np.sin(warp_theta) + x0, + # x + # ) + # y_warped = np.where( + # (x >= x0) & (x <= x1) & (np.abs(warp_radius) < np.inf), + # warp_radius * np.cos(warp_theta) + y * np.cos(warp_theta) - warp_radius, + # y + # ) + + ### Instead, use a simpler y "slide" approach + if len(x.shape) == len(x1.shape): + x1_rep = np.repeat(x1, x.shape[1], axis=1) # Usual solve + else: + x1_rep = x1[:,:,None] # Full field solve + x_mod = np.where( + x > x1_rep, + x1_rep, + x + ) + y_warped = np.where( + (np.abs(warp_radius) != np.inf ) & (x >= x0), + -warp_radius + np.sign(warp_radius) * np.sqrt(warp_radius**2 - (x_mod - x0)**2) + y_warped, + y_warped + ) # Return warped grid return x, y_warped, z From 845d3da36df070cbefd30a817295af087d31e6eb Mon Sep 17 00:00:00 2001 From: misi9170 Date: Thu, 23 May 2024 09:36:44 -0600 Subject: [PATCH 08/37] Enable other full_flow solvers. --- floris/core/solver.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/floris/core/solver.py b/floris/core/solver.py index 42479f734..a129ce8cf 100644 --- a/floris/core/solver.py +++ b/floris/core/solver.py @@ -724,6 +724,7 @@ def full_flow_cc_solver( turbine_coordinates=turbine_grid_farm.coordinates, turbine_diameters=turbine_grid_farm.rotor_diameters, wind_directions=turbine_grid_flow_field.wind_directions, + heterogeneous_inflow_config=turbine_grid_flow_field.heterogeneous_inflow_config, grid_resolution=3, ) turbine_grid_farm.expand_farm_properties( @@ -1384,6 +1385,7 @@ def full_flow_empirical_gauss_solver( turbine_coordinates=turbine_grid_farm.coordinates, turbine_diameters=turbine_grid_farm.rotor_diameters, wind_directions=turbine_grid_flow_field.wind_directions, + heterogeneous_inflow_config=turbine_grid_flow_field.heterogeneous_inflow_config, grid_resolution=3, ) turbine_grid_farm.expand_farm_properties( From 89d771da3d926fc25c898dc2b46b20dc16ca4524 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Wed, 26 Jun 2024 14:07:47 -0600 Subject: [PATCH 09/37] Update file path to run. --- examples/examples_heterogeneous/bulk_het.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/examples_heterogeneous/bulk_het.py b/examples/examples_heterogeneous/bulk_het.py index f2080c76a..5e225915c 100644 --- a/examples/examples_heterogeneous/bulk_het.py +++ b/examples/examples_heterogeneous/bulk_het.py @@ -6,7 +6,7 @@ # Get a test fi (single turbine at 0,0) -fmodel = FlorisModel("inputs/gch.yaml") +fmodel = FlorisModel("../inputs/gch.yaml") fmodel.set(layout_x=[0], layout_y=[0]) # Directly downstream at 270 degrees From 63fe52adda6c0f989a1bae00e1666e2b1c39cd72 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Wed, 26 Jun 2024 14:15:22 -0600 Subject: [PATCH 10/37] Ruff and isort. --- examples/examples_heterogeneous/bulk_het.py | 7 ++++--- examples/examples_heterogeneous/bulk_het_2.py | 2 +- floris/utilities.py | 5 +++-- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/examples/examples_heterogeneous/bulk_het.py b/examples/examples_heterogeneous/bulk_het.py index 5e225915c..0b5989d4e 100644 --- a/examples/examples_heterogeneous/bulk_het.py +++ b/examples/examples_heterogeneous/bulk_het.py @@ -1,8 +1,9 @@ -import numpy as np +from pathlib import Path + import matplotlib.pyplot as plt -from floris import FlorisModel +import numpy as np -from pathlib import Path +from floris import FlorisModel # Get a test fi (single turbine at 0,0) diff --git a/examples/examples_heterogeneous/bulk_het_2.py b/examples/examples_heterogeneous/bulk_het_2.py index 645e3b9cd..da1f49f1d 100644 --- a/examples/examples_heterogeneous/bulk_het_2.py +++ b/examples/examples_heterogeneous/bulk_het_2.py @@ -1,5 +1,5 @@ -import numpy as np import matplotlib.pyplot as plt +import numpy as np import floris.layout_visualization as layoutviz from floris import FlorisModel diff --git a/floris/utilities.py b/floris/utilities.py index 6bc88587c..1f2ce67d3 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -233,7 +233,7 @@ def warp_grid_for_wind_direction_heterogeneity( ): y_warped = y.copy() - + for i in range(wd_het_values.shape[1]-1): x0 = wd_het_x_points[:,i:i+1] @@ -271,7 +271,8 @@ def warp_grid_for_wind_direction_heterogeneity( ) y_warped = np.where( (np.abs(warp_radius) != np.inf ) & (x >= x0), - -warp_radius + np.sign(warp_radius) * np.sqrt(warp_radius**2 - (x_mod - x0)**2) + y_warped, + (-warp_radius + np.sign(warp_radius) * np.sqrt(warp_radius**2 - (x_mod - x0)**2) + + y_warped), y_warped ) From 46e28aa0d96afe908e338307fb28e5c5d4008bfe Mon Sep 17 00:00:00 2001 From: misi9170 Date: Thu, 25 Jul 2024 10:52:30 -0600 Subject: [PATCH 11/37] Initial infrastucture for turbine-specific layouts. Currently only in sequential_solver. --- floris/core/grid.py | 2 ++ floris/core/solver.py | 10 +++++++++- floris/utilities.py | 22 ++++++++++++++++++++++ 3 files changed, 33 insertions(+), 1 deletion(-) diff --git a/floris/core/grid.py b/floris/core/grid.py index fa01c600e..b11e8cbac 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -65,6 +65,8 @@ class Grid(ABC, BaseClass): z_sorted_inertial_frame: NDArrayFloat = field(init=False) cubature_weights: NDArrayFloat = field(init=False, default=None) + use_turbine_specific_layouts: NDArrayFloat = field(init=False, default=False) + @turbine_coordinates.validator def check_coordinates(self, instance: attrs.Attribute, value: np.ndarray) -> None: """ diff --git a/floris/core/solver.py b/floris/core/solver.py index a129ce8cf..c99b31ad7 100644 --- a/floris/core/solver.py +++ b/floris/core/solver.py @@ -25,7 +25,7 @@ ) from floris.core.wake_velocity.empirical_gauss import awc_added_wake_mixing from floris.type_dec import NDArrayFloat -from floris.utilities import cosd +from floris.utilities import cosd, update_model_args def calculate_area_overlap(wake_velocities, freestream_velocities, y_ngrid, z_ngrid): @@ -89,6 +89,14 @@ def sequential_solver( u_i = flow_field.u_sorted[:, i:i+1] v_i = flow_field.v_sorted[:, i:i+1] + if grid.use_turbine_specific_layouts: + deficit_model_args = update_model_args( + deficit_model_args, + grid.x_sorted_per_turbine[:, :, i], + grid.y_sorted_per_turbine[:, :, i], + grid.z_sorted_per_turbine[:, :, i], + ) + ct_i = thrust_coefficient( velocities=flow_field.u_sorted, turbulence_intensities=flow_field.turbulence_intensity_field_sorted, diff --git a/floris/utilities.py b/floris/utilities.py index 1f2ce67d3..d128565ab 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -393,3 +393,25 @@ def print_nested_dict(dictionary: Dict[str, Any], indent: int = 0) -> None: print_nested_dict(value, indent + 4) else: print(" " * (indent + 4) + str(value)) + +def update_model_args( + model_args: Dict[str, Any], + x: np.ndarray, + y: np.ndarray, + z: np.ndarray + ): + """ + Update the x_sorted, y_sorted, and z_sorted entries of the model_args dict. + """ + # Check that new x, y, z are the same length as the old + if (x.shape != model_args["x"].shape or + y.shape != model_args["y"].shape or + z.shape != model_args["z"].shape): + raise ValueError("Size of new x,y,z arrays must match existing ones.") + + # Update the x, y, z values + model_args["x"] = x + model_args["y"] = y + model_args["z"] = z + + return model_args From 4bd5e82438137c2277712019f6e27defbcc277d3 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Thu, 25 Jul 2024 11:04:50 -0600 Subject: [PATCH 12/37] Mock up with turbine locations simply repeated. --- floris/core/grid.py | 34 ++++++++++++---------------------- floris/core/solver.py | 6 +++--- 2 files changed, 15 insertions(+), 25 deletions(-) diff --git a/floris/core/grid.py b/floris/core/grid.py index b11e8cbac..cedbd4a87 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -66,6 +66,9 @@ class Grid(ABC, BaseClass): cubature_weights: NDArrayFloat = field(init=False, default=None) use_turbine_specific_layouts: NDArrayFloat = field(init=False, default=False) + x_sorted_per_turbine: NDArrayFloat = field(init=False) + y_sorted_per_turbine: NDArrayFloat = field(init=False) + z_sorted_per_turbine: NDArrayFloat = field(init=False) @turbine_coordinates.validator def check_coordinates(self, instance: attrs.Attribute, value: np.ndarray) -> None: @@ -191,28 +194,6 @@ def set_grid(self) -> None: self.turbine_coordinates, ) - if (self.heterogeneous_inflow_config is not None - and "bulk_wd_change" in self.heterogeneous_inflow_config): - wd_het_x_points = np.array(self.heterogeneous_inflow_config["bulk_wd_x"]) - wd_het_values = np.array(self.heterogeneous_inflow_config["bulk_wd_change"]) - # TODO: does the below work as expected? Do I need the y values too? - # coordinates_to_rotate = np.concatenate( - # ( - # wd_het_x_points[:,:,None], - # np.zeros((wd_het_x_points.shape[0], wd_het_x_points.shape[1], 2)) - # ), - # axis=2 - # )[None,:,:,:] - # wd_het_x_points = rotate_coordinates_rel_west( - # np.repeat(self.wind_directions[None,:], wd_het_x_points.shape[1], axis=0), - # coordinates_to_rotate, - # self.x_center_of_rotation, - # self.y_center_of_rotation, - # )[0][:, :, 0].T - x, y, z = warp_grid_for_wind_direction_heterogeneity( - x, y, z, wd_het_x_points + x.min(), wd_het_values - ) - # - **rloc** (*float, optional): A value, from 0 to 1, that determines # the width/height of the grid of points on the rotor as a ratio of # the rotor radius. @@ -271,6 +252,15 @@ def set_grid(self) -> None: self.y_sorted = np.take_along_axis(_y, self.sorted_indices, axis=1) self.z_sorted = np.take_along_axis(_z, self.sorted_indices, axis=1) + # if (self.heterogeneous_inflow_config is not None + # and "wind_directions" in self.heterogeneous_inflow_config): + if True: + self.use_turbine_specific_layouts = True + # set up new x, y, z coordinates as a mockup + self.x_sorted_per_turbine = np.repeat(self.x_sorted[:,:,None,:,:], self.n_turbines, axis=2) + self.y_sorted_per_turbine = np.repeat(self.y_sorted[:,:,None,:,:], self.n_turbines, axis=2) + self.z_sorted_per_turbine = np.repeat(self.z_sorted[:,:,None,:,:], self.n_turbines, axis=2) + # Now calculate grid coordinates in original frame (from 270 deg perspective) self.x_sorted_inertial_frame, self.y_sorted_inertial_frame, self.z_sorted_inertial_frame = \ reverse_rotate_coordinates_rel_west( diff --git a/floris/core/solver.py b/floris/core/solver.py index c99b31ad7..dc8fa47cb 100644 --- a/floris/core/solver.py +++ b/floris/core/solver.py @@ -92,9 +92,9 @@ def sequential_solver( if grid.use_turbine_specific_layouts: deficit_model_args = update_model_args( deficit_model_args, - grid.x_sorted_per_turbine[:, :, i], - grid.y_sorted_per_turbine[:, :, i], - grid.z_sorted_per_turbine[:, :, i], + grid.x_sorted_per_turbine[:, :, i, :, :], + grid.y_sorted_per_turbine[:, :, i, :, :], + grid.z_sorted_per_turbine[:, :, i, :, :], ) ct_i = thrust_coefficient( From b067dee5edb9cac9421539a91e255debc31f1294 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Thu, 25 Jul 2024 13:52:46 -0600 Subject: [PATCH 13/37] Placeholder for the layout modifier. --- floris/core/flow_field.py | 9 +++------ floris/core/grid.py | 21 ++++++++++++++------- floris/core/solver.py | 12 +++++++++--- floris/utilities.py | 18 ++++++++++++++++++ 4 files changed, 44 insertions(+), 16 deletions(-) diff --git a/floris/core/flow_field.py b/floris/core/flow_field.py index d7ea19616..f9e3bab2a 100644 --- a/floris/core/flow_field.py +++ b/floris/core/flow_field.py @@ -105,10 +105,8 @@ def heterogeneous_config_validator(self, instance: attrs.Attribute, value: dict return # Check that the correct keys are supplied for the heterogeneous_inflow_config dict - speed_het_keys = ["speed_multipliers", "x", "y", "z"] - bulk_wd_het_keys = ["bulk_wd_change", "bulk_wd_x"] - for k in value.keys(): - if k not in speed_het_keys and k not in bulk_wd_het_keys: + for k in ["speed_multipliers", "x", "y"]: + if k not in value.keys(): raise ValueError( "heterogeneous_inflow_config must contain entries for 'speed_multipliers'," f"'x', and 'y', with 'z' optional. Missing '{k}'." @@ -132,8 +130,7 @@ def het_map_validator(self, instance: attrs.Attribute, value: list | None) -> No def __attrs_post_init__(self) -> None: - if self.heterogeneous_inflow_config is not None and self._speed_heterogeneity: - print(self._speed_heterogeneity) + if self._speed_heterogeneity: self.generate_heterogeneous_wind_map() diff --git a/floris/core/grid.py b/floris/core/grid.py index cedbd4a87..2687504de 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -19,6 +19,7 @@ reverse_rotate_coordinates_rel_west, rotate_coordinates_rel_west, warp_grid_for_wind_direction_heterogeneity, + apply_wind_direction_heterogeneity, ) @@ -252,15 +253,21 @@ def set_grid(self) -> None: self.y_sorted = np.take_along_axis(_y, self.sorted_indices, axis=1) self.z_sorted = np.take_along_axis(_z, self.sorted_indices, axis=1) - # if (self.heterogeneous_inflow_config is not None - # and "wind_directions" in self.heterogeneous_inflow_config): - if True: + if (self.heterogeneous_inflow_config is not None + and "wind_directions" in self.heterogeneous_inflow_config): self.use_turbine_specific_layouts = True # set up new x, y, z coordinates as a mockup - self.x_sorted_per_turbine = np.repeat(self.x_sorted[:,:,None,:,:], self.n_turbines, axis=2) - self.y_sorted_per_turbine = np.repeat(self.y_sorted[:,:,None,:,:], self.n_turbines, axis=2) - self.z_sorted_per_turbine = np.repeat(self.z_sorted[:,:,None,:,:], self.n_turbines, axis=2) - + self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \ + apply_wind_direction_heterogeneity( + self.x_sorted, + self.y_sorted, + self.z_sorted, + None, + None, + None, + None, + self.n_turbines + ) # Now calculate grid coordinates in original frame (from 270 deg perspective) self.x_sorted_inertial_frame, self.y_sorted_inertial_frame, self.z_sorted_inertial_frame = \ reverse_rotate_coordinates_rel_west( diff --git a/floris/core/solver.py b/floris/core/solver.py index dc8fa47cb..03047309d 100644 --- a/floris/core/solver.py +++ b/floris/core/solver.py @@ -92,9 +92,15 @@ def sequential_solver( if grid.use_turbine_specific_layouts: deficit_model_args = update_model_args( deficit_model_args, - grid.x_sorted_per_turbine[:, :, i, :, :], - grid.y_sorted_per_turbine[:, :, i, :, :], - grid.z_sorted_per_turbine[:, :, i, :, :], + grid.x_sorted_per_turbine[:, :, :, :, i], + grid.y_sorted_per_turbine[:, :, :, :, i], + grid.z_sorted_per_turbine[:, :, :, :, i], + ) + deflection_model_args = update_model_args( + deflection_model_args, + grid.x_sorted_per_turbine[:, :, :, :, i], + grid.y_sorted_per_turbine[:, :, :, :, i], + grid.z_sorted_per_turbine[:, :, :, :, i], ) ct_i = thrust_coefficient( diff --git a/floris/utilities.py b/floris/utilities.py index d128565ab..d5361b9ee 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -279,6 +279,24 @@ def warp_grid_for_wind_direction_heterogeneity( # Return warped grid return x, y_warped, z +def apply_wind_direction_heterogeneity( + x, + y, + z, + wind_direction_x_points, + wind_direction_y_points, + wind_direction_z_points, + wind_direction_values, + n_turbines + ): + + x_per_turbine = np.repeat(x[:,:,:,:,None], n_turbines, axis=4) + y_per_turbine = np.repeat(y[:,:,:,:,None], n_turbines, axis=4) + z_per_turbine = np.repeat(z[:,:,:,:,None], n_turbines, axis=4) + + # Return warped grid + return x_per_turbine, y_per_turbine, z_per_turbine + class Loader(yaml.SafeLoader): From 34e8b5393b62b67953bfefdbdee026142cd99fee Mon Sep 17 00:00:00 2001 From: misi9170 Date: Thu, 25 Jul 2024 14:37:03 -0600 Subject: [PATCH 14/37] Add basic example; wind_direcitons not used, simply hardcoded rearranging as proof of concept. --- examples/examples_heterogeneous/bulk_het_3.py | 33 +++++++++++++++++++ floris/utilities.py | 12 ++++++- 2 files changed, 44 insertions(+), 1 deletion(-) create mode 100644 examples/examples_heterogeneous/bulk_het_3.py diff --git a/examples/examples_heterogeneous/bulk_het_3.py b/examples/examples_heterogeneous/bulk_het_3.py new file mode 100644 index 000000000..61b4f3fd9 --- /dev/null +++ b/examples/examples_heterogeneous/bulk_het_3.py @@ -0,0 +1,33 @@ +import matplotlib.pyplot as plt +import numpy as np + +import floris.layout_visualization as layoutviz +from floris import FlorisModel +from floris.flow_visualization import visualize_cut_plane + + +# Get a test fi (single turbine at 0,0) +fmodel = FlorisModel("../inputs/gch.yaml") +fmodel.set(layout_x=[0, 0, 600, 600], layout_y=[0, -100, 0, -100]) + +fmodel.set( + wind_directions=[270.0], + wind_speeds=[8.0], + turbulence_intensities=[0.06], +) +fmodel.run() +P_wo_het = fmodel.get_turbine_powers()/1e6 + +fmodel.set( + heterogeneous_inflow_config={ + "wind_directions": [270.0, 280.0, 280.0, 290.0], + "x": [-1000.0, -1000.0, 1000.0, 1000.0], + "y": [-1000.0, 1000.0, -1000.0, 1000.0], + "speed_multipliers":np.array([[1.0, 1.0, 1.0, 1.0]])#np.array([[1.0, 1.2, 1.2, 1.4]]) + }, +) +fmodel.run() +P_w_het = fmodel.get_turbine_powers()/1e6 + +print("Difference (MW):", P_w_het - P_wo_het) + diff --git a/floris/utilities.py b/floris/utilities.py index d5361b9ee..e36835938 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -290,11 +290,21 @@ def apply_wind_direction_heterogeneity( n_turbines ): + # Initialize storage x_per_turbine = np.repeat(x[:,:,:,:,None], n_turbines, axis=4) y_per_turbine = np.repeat(y[:,:,:,:,None], n_turbines, axis=4) z_per_turbine = np.repeat(z[:,:,:,:,None], n_turbines, axis=4) - # Return warped grid + # Compute new relative locations + # TODO + # TEMP________ + # Place turbine 2 directly downstream of turbine 1, as seen by turbine 1 + y_per_turbine[0,2,:,:,1] = y_per_turbine[0,1,:,:,1] + # Place turbine 3 out of the direct path of turbine 1, as seen by turbine 1 + y_per_turbine[0,3,:,:,1] = 2*y_per_turbine[0,3,:,:,1] + # END TEMP____ + + # Return full set of locations return x_per_turbine, y_per_turbine, z_per_turbine From 66ef9fed1193a67d082aa3a0f2847dd04e233df6 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Mon, 29 Jul 2024 12:34:20 -0600 Subject: [PATCH 15/37] Infrastructure now in place and running for full_flow_sequential_solver --- examples/examples_heterogeneous/bulk_het_3.py | 54 ++++++++++++++++--- floris/core/grid.py | 22 ++++++-- floris/core/solver.py | 14 +++++ floris/floris_model.py | 25 +++++---- floris/utilities.py | 17 ++++-- 5 files changed, 107 insertions(+), 25 deletions(-) diff --git a/examples/examples_heterogeneous/bulk_het_3.py b/examples/examples_heterogeneous/bulk_het_3.py index 61b4f3fd9..c9e00b80e 100644 --- a/examples/examples_heterogeneous/bulk_het_3.py +++ b/examples/examples_heterogeneous/bulk_het_3.py @@ -5,6 +5,7 @@ from floris import FlorisModel from floris.flow_visualization import visualize_cut_plane +visualize = True # Get a test fi (single turbine at 0,0) fmodel = FlorisModel("../inputs/gch.yaml") @@ -18,16 +19,53 @@ fmodel.run() P_wo_het = fmodel.get_turbine_powers()/1e6 -fmodel.set( - heterogeneous_inflow_config={ - "wind_directions": [270.0, 280.0, 280.0, 290.0], - "x": [-1000.0, -1000.0, 1000.0, 1000.0], - "y": [-1000.0, 1000.0, -1000.0, 1000.0], - "speed_multipliers":np.array([[1.0, 1.0, 1.0, 1.0]])#np.array([[1.0, 1.2, 1.2, 1.4]]) - }, -) +het_inflow_config = { + "wind_directions": [270.0, 280.0, 280.0, 290.0], + "x": np.array([-1000.0, -1000.0, 1000.0, 1000.0]), + "y": np.array([-1000.0, 1000.0, -1000.0, 1000.0]), + "speed_multipliers":np.array([[1.0, 1.0, 1.0, 1.0]])#np.array([[1.0, 1.2, 1.2, 1.4]]) +} + +fmodel.set(heterogeneous_inflow_config=het_inflow_config) fmodel.run() P_w_het = fmodel.get_turbine_powers()/1e6 print("Difference (MW):", P_w_het - P_wo_het) +if visualize: + fig, ax = plt.subplots(2, 1) + fmodel.set( + heterogeneous_inflow_config={ + #"wind_directions": [270.0, 270.0, 270.0, 270.0], + "x": np.array([-1000.0, -1000.0, 1000.0, 1000.0]), + "y": np.array([-1000.0, 1000.0, -1000.0, 1000.0]), + "speed_multipliers":np.array([[1.0, 1.0, 1.0, 1.0]])#np.array([[1.0, 1.2, 1.2, 1.4]]) + } + ) + fmodel.run() + horizontal_plane = fmodel.calculate_horizontal_plane( + x_resolution=200, + y_resolution=100, + height=90.0, + ) + visualize_cut_plane( + horizontal_plane, + ax=ax[0], + label_contours=False, + title="Without WD het" + ) + + fmodel.set(heterogeneous_inflow_config=het_inflow_config) + horizontal_plane = fmodel.calculate_horizontal_plane( + x_resolution=200, + y_resolution=100, + height=90.0, + ) + visualize_cut_plane( + horizontal_plane, + ax=ax[1], + label_contours=False, + title="With WD het" + ) + + plt.show() diff --git a/floris/core/grid.py b/floris/core/grid.py index 2687504de..c3fc5df43 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -67,9 +67,9 @@ class Grid(ABC, BaseClass): cubature_weights: NDArrayFloat = field(init=False, default=None) use_turbine_specific_layouts: NDArrayFloat = field(init=False, default=False) - x_sorted_per_turbine: NDArrayFloat = field(init=False) - y_sorted_per_turbine: NDArrayFloat = field(init=False) - z_sorted_per_turbine: NDArrayFloat = field(init=False) + x_sorted_per_turbine: NDArrayFloat = field(init=False, default=None) + y_sorted_per_turbine: NDArrayFloat = field(init=False, default=None) + z_sorted_per_turbine: NDArrayFloat = field(init=False, default=None) @turbine_coordinates.validator def check_coordinates(self, instance: attrs.Attribute, value: np.ndarray) -> None: @@ -268,6 +268,7 @@ def set_grid(self) -> None: None, self.n_turbines ) + # Now calculate grid coordinates in original frame (from 270 deg perspective) self.x_sorted_inertial_frame, self.y_sorted_inertial_frame, self.z_sorted_inertial_frame = \ reverse_rotate_coordinates_rel_west( @@ -627,6 +628,21 @@ def set_grid(self) -> None: self.y_sorted = y_points[None, :, :, :] self.z_sorted = z_points[None, :, :, :] + if (self.heterogeneous_inflow_config is not None + and "wind_directions" in self.heterogeneous_inflow_config): + self.use_turbine_specific_layouts = True + # set up new x, y, z coordinates as a mockup + self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \ + apply_wind_direction_heterogeneity( + self.x_sorted, + self.y_sorted, + self.z_sorted, + None, + None, + None, + None, + self.n_turbines + ) # Now calculate grid coordinates in original frame (from 270 deg perspective) self.x_sorted_inertial_frame, self.y_sorted_inertial_frame, self.z_sorted_inertial_frame = \ diff --git a/floris/core/solver.py b/floris/core/solver.py index 03047309d..8668e5c4f 100644 --- a/floris/core/solver.py +++ b/floris/core/solver.py @@ -342,6 +342,20 @@ def full_flow_sequential_solver( u_i = turbine_grid_flow_field.u_sorted[:, i:i+1] v_i = turbine_grid_flow_field.v_sorted[:, i:i+1] + if flow_field_grid.use_turbine_specific_layouts: + deficit_model_args = update_model_args( + deficit_model_args, + flow_field_grid.x_sorted_per_turbine[:, :, :, :, i], + flow_field_grid.y_sorted_per_turbine[:, :, :, :, i], + flow_field_grid.z_sorted_per_turbine[:, :, :, :, i], + ) + deflection_model_args = update_model_args( + deflection_model_args, + flow_field_grid.x_sorted_per_turbine[:, :, :, :, i], + flow_field_grid.y_sorted_per_turbine[:, :, :, :, i], + flow_field_grid.z_sorted_per_turbine[:, :, :, :, i], + ) + ct_i = thrust_coefficient( velocities=turbine_grid_flow_field.u_sorted, turbulence_intensities=turbine_grid_flow_field.turbulence_intensity_field_sorted, diff --git a/floris/floris_model.py b/floris/floris_model.py index 73d1f1fd5..de84d7d30 100644 --- a/floris/floris_model.py +++ b/floris/floris_model.py @@ -969,18 +969,23 @@ def set_for_viz(self, findex: int, solver_settings: dict) -> None: """ # If not None, set the heterogeneous inflow configuration + # TODO: add wind direction info here if self.core.flow_field.heterogeneous_inflow_config is not None: - speed_het_keys = ["speed_multipliers", "x", "y", "z"] - bulk_wd_het_keys = ["bulk_wd_change", "bulk_wd_x"] - heterogeneous_inflow_config = {} - for k in self.core.flow_field.heterogeneous_inflow_config.keys(): - if ( - (k in speed_het_keys - or k in bulk_wd_het_keys) - and self.core.flow_field.heterogeneous_inflow_config[k] is not None - ): + heterogeneous_inflow_config = { + 'x': self.core.flow_field.heterogeneous_inflow_config['x'], + 'y': self.core.flow_field.heterogeneous_inflow_config['y'], + 'speed_multipliers': + self.core.flow_field.heterogeneous_inflow_config['speed_multipliers'][findex:findex+1], + } + for k in ["z"]: # Optional, independent of findex + if k in self.core.flow_field.heterogeneous_inflow_config: heterogeneous_inflow_config[k] = ( - self.core.flow_field.heterogeneous_inflow_config[k][findex:findex+1] + self.core.flow_field.heterogeneous_inflow_config[k] + ) + for k in ["wind_directions"]: # Optional, dependent on findex + if k in self.core.flow_field.heterogeneous_inflow_config: + heterogeneous_inflow_config[k] = ( + self.core.flow_field.heterogeneous_inflow_config[k][findex:findex+1], ) else: heterogeneous_inflow_config = None diff --git a/floris/utilities.py b/floris/utilities.py index e36835938..81f4fc10b 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -298,10 +298,19 @@ def apply_wind_direction_heterogeneity( # Compute new relative locations # TODO # TEMP________ - # Place turbine 2 directly downstream of turbine 1, as seen by turbine 1 - y_per_turbine[0,2,:,:,1] = y_per_turbine[0,1,:,:,1] - # Place turbine 3 out of the direct path of turbine 1, as seen by turbine 1 - y_per_turbine[0,3,:,:,1] = 2*y_per_turbine[0,3,:,:,1] + option = 3 + if option == 1: + # Place turbine 2 directly downstream of turbine 1, as seen by turbine 1 + y_per_turbine[0,2,:,:,1] = y_per_turbine[0,1,:,:,1] + # Place turbine 3 out of the direct path of turbine 1, as seen by turbine 1 + y_per_turbine[0,3,:,:,1] = 2*y_per_turbine[0,3,:,:,1] + elif option == 2: + #y_per_turbine[0,:,:,:,0] = y_per_turbine[0,:,:,:,0] - 0.2*x_per_turbine[0,:,:,:,0] + #y_per_turbine[0,:,:,:,1] = y_per_turbine[0,:,:,:,1] + 0.2*x_per_turbine[0,:,:,:,1] + y_per_turbine[0,:,:,:,0] = y_per_turbine[0,:,:,:,0] + 0.2*x_per_turbine[0,:,:,:,0] + else: + pass # make no change + # END TEMP____ # Return full set of locations From 3d73fe39c33c1bf7eb401b4ce898c8437b7be59e Mon Sep 17 00:00:00 2001 From: misi9170 Date: Mon, 29 Jul 2024 16:46:14 -0600 Subject: [PATCH 16/37] Working through a 2D example case. --- examples/examples_heterogeneous/bulk_het_4.py | 76 +++++++++++ floris/core/grid.py | 23 ++-- floris/utilities.py | 118 +++++++++++++++++- 3 files changed, 205 insertions(+), 12 deletions(-) create mode 100644 examples/examples_heterogeneous/bulk_het_4.py diff --git a/examples/examples_heterogeneous/bulk_het_4.py b/examples/examples_heterogeneous/bulk_het_4.py new file mode 100644 index 000000000..748e8546c --- /dev/null +++ b/examples/examples_heterogeneous/bulk_het_4.py @@ -0,0 +1,76 @@ +import matplotlib.pyplot as plt +import numpy as np + +import floris.layout_visualization as layoutviz +from floris import FlorisModel +from floris.flow_visualization import visualize_cut_plane + +visualize = True + +# Get a test fi (single turbine at 0,0) +fmodel = FlorisModel("../inputs/gch.yaml") +fmodel.set(layout_x=[0, 0, 600, 600], layout_y=[0, -100, 0, -100]) + +fmodel.set( + wind_directions=[270.0], + wind_speeds=[8.0], + turbulence_intensities=[0.06], + wind_shear=0.0 +) +fmodel.run() +P_wo_het = fmodel.get_turbine_powers()/1e6 + +het_inflow_config = { + "wind_directions": [270.0, 280.0, 280.0, 290.0], # To trigger the needed code + "x": np.array([-1000.0, -1000.0, 1000.0, 1000.0]), + "y": np.array([-1000.0, 1000.0, -1000.0, 1000.0]), + #"z": np.array([90.0, 90.0, 90.0, 91.0]), + "speed_multipliers": np.array([[1.0, 1.0, 1.0, 1.0]]), + "u": np.array([[8.0, 8.0, 8.0, 8.0]]), + "v": np.array([[-2.0, 0.0, -2.0, 0.0]]), + #"w": np.array([[0.0, 0.0, 0.0, 0.0]]), +} + +fmodel.set(heterogeneous_inflow_config=het_inflow_config) +fmodel.run() +P_w_het = fmodel.get_turbine_powers()/1e6 + +print("Difference (MW):", P_w_het - P_wo_het) + +if visualize: + fig, ax = plt.subplots(2, 1) + fmodel.set( + heterogeneous_inflow_config={ + #"wind_directions": [270.0, 270.0, 270.0, 270.0], + "x": np.array([-1000.0, -1000.0, 1000.0, 1000.0]), + "y": np.array([-1000.0, 1000.0, -1000.0, 1000.0]), + "speed_multipliers":np.array([[1.0, 1.0, 1.0, 1.0]])#np.array([[1.0, 1.2, 1.2, 1.4]]) + } + ) + fmodel.run() + horizontal_plane = fmodel.calculate_horizontal_plane( + x_resolution=200, + y_resolution=100, + height=90.0, + ) + visualize_cut_plane( + horizontal_plane, + ax=ax[0], + label_contours=False, + title="Without WD het" + ) + + fmodel.set(heterogeneous_inflow_config=het_inflow_config) + horizontal_plane = fmodel.calculate_horizontal_plane( + x_resolution=200, + y_resolution=100, + height=90.0, + ) + visualize_cut_plane( + horizontal_plane, + ax=ax[1], + label_contours=False, + title="With WD het" + ) + + plt.show() diff --git a/floris/core/grid.py b/floris/core/grid.py index c3fc5df43..d3a4ba505 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -262,10 +262,12 @@ def set_grid(self) -> None: self.x_sorted, self.y_sorted, self.z_sorted, - None, - None, - None, - None, + self.heterogeneous_inflow_config["x"], + self.heterogeneous_inflow_config["y"], + self.heterogeneous_inflow_config["z"] if "z" in self.heterogeneous_inflow_config else None, + self.heterogeneous_inflow_config["u"], + self.heterogeneous_inflow_config["v"], + self.heterogeneous_inflow_config["w"] if "w" in self.heterogeneous_inflow_config else None, self.n_turbines ) @@ -637,13 +639,14 @@ def set_grid(self) -> None: self.x_sorted, self.y_sorted, self.z_sorted, - None, - None, - None, - None, + self.heterogeneous_inflow_config["x"], + self.heterogeneous_inflow_config["y"], + self.heterogeneous_inflow_config["z"], + self.heterogeneous_inflow_config["u"], + self.heterogeneous_inflow_config["v"], + self.heterogeneous_inflow_config["w"], self.n_turbines ) - # Now calculate grid coordinates in original frame (from 270 deg perspective) self.x_sorted_inertial_frame, self.y_sorted_inertial_frame, self.z_sorted_inertial_frame = \ reverse_rotate_coordinates_rel_west( @@ -673,7 +676,7 @@ class PointsGrid(Grid): Args: turbine_coordinates (:py:obj:`NDArrayFloat`): Not used for PointsGrid, but required for the `Grid` super-class. - turbine_diameters (:py:obj:`NDArrayFloat`): Not used for PointsGrid, but + turbine_diameters (:py:obj:`NDArrayFloat`): Not used for PointsGrid, but required for the `Grid` super-class. wind_directions (:py:obj:`NDArrayFloat`): Wind directions supplied by the user. grid_resolution (:py:obj:`int` | :py:obj:`Iterable(int,)`): Not used for PointsGrid, but diff --git a/floris/utilities.py b/floris/utilities.py index 81f4fc10b..8b923b889 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -14,6 +14,8 @@ import numpy as np import yaml from attrs import define, field +from scipy.interpolate import LinearNDInterpolator +from scipy.integrate import odeint from floris.type_dec import floris_array_converter, NDArrayFloat @@ -286,17 +288,20 @@ def apply_wind_direction_heterogeneity( wind_direction_x_points, wind_direction_y_points, wind_direction_z_points, - wind_direction_values, + wind_direction_u_values, + wind_direction_v_values, + wind_direction_w_values, n_turbines ): + compute_streamline = compute_streamline_2D if wind_direction_z_points is None else compute_streamline_3D + # Initialize storage x_per_turbine = np.repeat(x[:,:,:,:,None], n_turbines, axis=4) y_per_turbine = np.repeat(y[:,:,:,:,None], n_turbines, axis=4) z_per_turbine = np.repeat(z[:,:,:,:,None], n_turbines, axis=4) # Compute new relative locations - # TODO # TEMP________ option = 3 if option == 1: @@ -308,6 +313,25 @@ def apply_wind_direction_heterogeneity( #y_per_turbine[0,:,:,:,0] = y_per_turbine[0,:,:,:,0] - 0.2*x_per_turbine[0,:,:,:,0] #y_per_turbine[0,:,:,:,1] = y_per_turbine[0,:,:,:,1] + 0.2*x_per_turbine[0,:,:,:,1] y_per_turbine[0,:,:,:,0] = y_per_turbine[0,:,:,:,0] + 0.2*x_per_turbine[0,:,:,:,0] + elif option == 3: + # TODO: compute streamlines; convert to point movements. + # TODO: do I need an outer loop over the findices? Maybe; ideally not. + for tindex in range(n_turbines): + # 1. Compute the streamline + streamline = compute_streamline( + wind_direction_x_points, + wind_direction_y_points, + wind_direction_z_points, + wind_direction_u_values, + wind_direction_v_values, + wind_direction_w_values, + np.array([-100.0, 50.0]) # Placeholder + ) + import ipdb; ipdb.set_trace() + # 2. Compute the point movements + # 3. Apply the point movements to x_per_turbine, y_per_turbine, z_per_turbine + y_per_turbine[:,:,:,tindex] = 0 # Placeholder + else: pass # make no change @@ -316,6 +340,96 @@ def apply_wind_direction_heterogeneity( # Return full set of locations return x_per_turbine, y_per_turbine, z_per_turbine +def compute_streamline_3D(x, y, z, u, v, w, xyz_0): + """ + Compute streamline starting at xyz_0. + """ + s = np.linspace(0, 1, 1000) # parametric variable + + scale = np.array([x.max() - x.min(), y.max() - y.min(), z.max() - z.min()]) + offset = np.array([(x.max() + x.min())/2, (y.max() + y.min())/2, (z.max() + z.min())/2]) + # Need to compute the offset too + + # Collect + xyz = np.array([x, y, z]).T + uvw = np.concatenate((u, v, w), axis=0).T + + # Normalize + xyz = (xyz - offset) / scale + xyz_0 = (xyz_0 - offset) / scale + # TODO: Loop over findices! Maybe not necessary---can I compute all streamlines at once?? + # Would it help? Currently assuming a single findex. + interp = LinearNDInterpolator(xyz, uvw, fill_value=0.0) + + def velocity_field_interpolate_reg(xyz_, s): + # if (xyz_ <= -0.5).any() or (xyz_ >= 0.5).any(): + # print("hmm") + # return np.array([0.0, 0.0, 0.0]) + # else: + # return interp(xyz_)[0] + return interp(xyz_)[0] + + streamline = odeint(velocity_field_interpolate_reg, xyz_0, s) # 1000 x 3 + + # Denormalize + streamline = streamline * scale + offset + + return streamline + +def compute_streamline_2D(x, y, z, u, v, w, xy_0): + """ + Compute streamline starting at xy_0. + z and w will be ignored. + """ + s = np.linspace(0, 1, 1000) # parametric variable + + scale = np.array([x.max() - x.min(), y.max() - y.min()]) + offset = np.array([(x.max() + x.min())/2, (y.max() + y.min())/2]) + # Need to compute the offset too + + # Collect + xy = np.array([x, y]).T + uv = np.concatenate((u, v), axis=0).T + + # Normalize + xy = (xy - offset) / scale + xy_0 = (xy_0 - offset) / scale + # TODO: Loop over findices! Maybe not necessary---can I compute all streamlines at once?? + # Would it help? Currently assuming a single findex. + + interp = LinearNDInterpolator(xy, uv) + + def velocity_field_interpolate_reg(xy_, s): + if (xy_ <= -0.5).any() or (xy_ >= 0.5).any(): + return np.array([0, 0]) + else: + return interp(xy_)[0] + #return interp(xy_)[0] + + streamline = odeint(velocity_field_interpolate_reg, xy_0, s) # 1000 x 2 + + # Determine point of exit from specified domain + if (streamline < -0.5).any(): + low_exit = np.argwhere(streamline < -0.5)[0,0] + else: + low_exit = 1000 + if (streamline > 0.5).any(): + high_exit = np.argwhere(streamline > 0.5)[0,0] + else: + high_exit = 1000 + exit_index = np.minimum(low_exit, high_exit) + if exit_index == 0: + raise ValueError("Streamline could not be propagated inside the domain.") + if exit_index == 1000: + raise ValueError("Streamline did not reach the edge of the domain.") + + # Copy points after exit + streamline[exit_index:] = streamline[exit_index-1] + + # Denormalize + streamline = streamline * scale + offset + + return streamline class Loader(yaml.SafeLoader): From 9746ec821703def73b9c05aa593003283fca6fd6 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Wed, 31 Jul 2024 13:31:29 -0600 Subject: [PATCH 17/37] Viz by sweep working for single findex case. --- examples/examples_heterogeneous/bulk_het_4.py | 3 +- .../examples_heterogeneous/viz_by_sweep.py | 58 +++++++++++++++++ floris/core/grid.py | 20 +++--- floris/floris_model.py | 5 +- floris/utilities.py | 62 +++++++++++++++---- 5 files changed, 123 insertions(+), 25 deletions(-) create mode 100644 examples/examples_heterogeneous/viz_by_sweep.py diff --git a/examples/examples_heterogeneous/bulk_het_4.py b/examples/examples_heterogeneous/bulk_het_4.py index 748e8546c..326b9e33f 100644 --- a/examples/examples_heterogeneous/bulk_het_4.py +++ b/examples/examples_heterogeneous/bulk_het_4.py @@ -21,14 +21,13 @@ P_wo_het = fmodel.get_turbine_powers()/1e6 het_inflow_config = { - "wind_directions": [270.0, 280.0, 280.0, 290.0], # To trigger the needed code "x": np.array([-1000.0, -1000.0, 1000.0, 1000.0]), "y": np.array([-1000.0, 1000.0, -1000.0, 1000.0]), #"z": np.array([90.0, 90.0, 90.0, 91.0]), "speed_multipliers": np.array([[1.0, 1.0, 1.0, 1.0]]), "u": np.array([[8.0, 8.0, 8.0, 8.0]]), "v": np.array([[-2.0, 0.0, -2.0, 0.0]]), - #"w": np.array([[0.0, 0.0, 0.0, 0.0]]), + #"v": np.array([[0.0, 0.0, 0.0, 0.0]]), } fmodel.set(heterogeneous_inflow_config=het_inflow_config) diff --git a/examples/examples_heterogeneous/viz_by_sweep.py b/examples/examples_heterogeneous/viz_by_sweep.py new file mode 100644 index 000000000..5f846871a --- /dev/null +++ b/examples/examples_heterogeneous/viz_by_sweep.py @@ -0,0 +1,58 @@ +"""Example: Visualize flow by sweeping turbines + +Demonstrate the use calculate_horizontal_plane_with_turbines + +""" + +import matplotlib.pyplot as plt +import numpy as np + +import floris.flow_visualization as flowviz +from floris import FlorisModel + + +fmodel = FlorisModel("../inputs/gch.yaml") + +# # Some wake models may not yet have a visualization method included, for these cases can use +# # a slower version which scans a turbine model to produce the horizontal flow +x = np.array([-100, 0, 100, 200, 300, 400, 500, 600]) +v = -np.array([[0, 1, 2, 3, 4, 5, 6, 7]]) + +het_inflow_config = { + "x": np.repeat(x, 2), + "y": np.tile(np.array([-100.0, 100.0]), 8), + #"z": np.array([90.0, 90.0, 90.0, 91.0]), + "speed_multipliers": 1.0*np.ones((1, 16)), + "u": 8.0*np.ones((1, v.shape[1]*2)), + "v": np.repeat(v, 2, axis=1) + #"v": np.array([[0.0, 0.0, 0.0, 0.0]]), +} + +# Set a 2 turbine layout +fmodel.set( + layout_x=[0], + layout_y=[0], + wind_directions=[270], + wind_speeds=[8], + turbulence_intensities=[0.06], + heterogeneous_inflow_config=het_inflow_config, +) + +horizontal_plane_scan_turbine = flowviz.calculate_horizontal_plane_with_turbines( + fmodel, + x_resolution=20, + y_resolution=10, + x_bounds=(-100, 500), + y_bounds=(-100, 100), +) + +fig, ax = plt.subplots(figsize=(10, 4)) +flowviz.visualize_cut_plane( + horizontal_plane_scan_turbine, + ax=ax, + label_contours=True, + title="Horizontal (coarse turbine scan method)", +) + + +plt.show() diff --git a/floris/core/grid.py b/floris/core/grid.py index d3a4ba505..6f98add39 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -254,7 +254,7 @@ def set_grid(self) -> None: self.z_sorted = np.take_along_axis(_z, self.sorted_indices, axis=1) if (self.heterogeneous_inflow_config is not None - and "wind_directions" in self.heterogeneous_inflow_config): + and "u" in self.heterogeneous_inflow_config): self.use_turbine_specific_layouts = True # set up new x, y, z coordinates as a mockup self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \ @@ -262,12 +262,12 @@ def set_grid(self) -> None: self.x_sorted, self.y_sorted, self.z_sorted, - self.heterogeneous_inflow_config["x"], - self.heterogeneous_inflow_config["y"], - self.heterogeneous_inflow_config["z"] if "z" in self.heterogeneous_inflow_config else None, - self.heterogeneous_inflow_config["u"], - self.heterogeneous_inflow_config["v"], - self.heterogeneous_inflow_config["w"] if "w" in self.heterogeneous_inflow_config else None, + np.array(self.heterogeneous_inflow_config["x"]), + np.array(self.heterogeneous_inflow_config["y"]), + np.array(self.heterogeneous_inflow_config["z"]) if "z" in self.heterogeneous_inflow_config else None, + np.array(self.heterogeneous_inflow_config["u"]), + np.array(self.heterogeneous_inflow_config["v"]), + np.array(self.heterogeneous_inflow_config["w"]) if "w" in self.heterogeneous_inflow_config else None, self.n_turbines ) @@ -631,7 +631,7 @@ def set_grid(self) -> None: self.z_sorted = z_points[None, :, :, :] if (self.heterogeneous_inflow_config is not None - and "wind_directions" in self.heterogeneous_inflow_config): + and "u" in self.heterogeneous_inflow_config): self.use_turbine_specific_layouts = True # set up new x, y, z coordinates as a mockup self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \ @@ -641,10 +641,10 @@ def set_grid(self) -> None: self.z_sorted, self.heterogeneous_inflow_config["x"], self.heterogeneous_inflow_config["y"], - self.heterogeneous_inflow_config["z"], + self.heterogeneous_inflow_config["z"] if "z" in self.heterogeneous_inflow_config else None, self.heterogeneous_inflow_config["u"], self.heterogeneous_inflow_config["v"], - self.heterogeneous_inflow_config["w"], + self.heterogeneous_inflow_config["w"] if "w" in self.heterogeneous_inflow_config else None, self.n_turbines ) # Now calculate grid coordinates in original frame (from 270 deg perspective) diff --git a/floris/floris_model.py b/floris/floris_model.py index de84d7d30..3cf8ff8df 100644 --- a/floris/floris_model.py +++ b/floris/floris_model.py @@ -982,10 +982,11 @@ def set_for_viz(self, findex: int, solver_settings: dict) -> None: heterogeneous_inflow_config[k] = ( self.core.flow_field.heterogeneous_inflow_config[k] ) - for k in ["wind_directions"]: # Optional, dependent on findex + for k in ["wind_directions", "u", "v", "w"]: # Optional, dependent on findex + # TODO: add findex dimension back in! if k in self.core.flow_field.heterogeneous_inflow_config: heterogeneous_inflow_config[k] = ( - self.core.flow_field.heterogeneous_inflow_config[k][findex:findex+1], + self.core.flow_field.heterogeneous_inflow_config[k]#[findex:findex+1], ) else: heterogeneous_inflow_config = None diff --git a/floris/utilities.py b/floris/utilities.py index 8b923b889..cb3555b21 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -294,7 +294,8 @@ def apply_wind_direction_heterogeneity( n_turbines ): - compute_streamline = compute_streamline_2D if wind_direction_z_points is None else compute_streamline_3D + compute_streamline = compute_streamline_2D# if wind_direction_z_points is None else compute_streamline_3D + shift_points_by_streamline = shift_points_by_streamline_2D# if wind_direction_z_points is None else shift_points_by_streamline_3D # Initialize storage x_per_turbine = np.repeat(x[:,:,:,:,None], n_turbines, axis=4) @@ -325,13 +326,16 @@ def apply_wind_direction_heterogeneity( wind_direction_u_values, wind_direction_v_values, wind_direction_w_values, - np.array([-100.0, 50.0]) # Placeholder + np.array([x[0,tindex,:,:].mean(), y[0,tindex,:,:].mean()]) ) - import ipdb; ipdb.set_trace() + # 2. Compute the point movements - # 3. Apply the point movements to x_per_turbine, y_per_turbine, z_per_turbine - y_per_turbine[:,:,:,tindex] = 0 # Placeholder - + x_shifted, y_shifted, z_shifted = shift_points_by_streamline(streamline, x, y, z) + + # 3. Apply to per_turbine values + x_per_turbine[:,:,:,:,tindex] = x_shifted + y_per_turbine[:,:,:,:,tindex] = y_shifted + z_per_turbine[:,:,:,:,tindex] = z_shifted else: pass # make no change @@ -376,12 +380,12 @@ def velocity_field_interpolate_reg(xyz_, s): return streamline -def compute_streamline_2D(x, y, z, u, v, w, xy_0): +def compute_streamline_2D(x, y, z, u, v, w, xy_0, n_points=1000): """ Compute streamline starting at xy_0. z and w will be ignored. """ - s = np.linspace(0, 1, 1000) # parametric variable + s = np.linspace(0, 1, n_points) # parametric variable scale = np.array([x.max() - x.min(), y.max() - y.min()]) offset = np.array([(x.max() + x.min())/2, (y.max() + y.min())/2]) @@ -400,13 +404,12 @@ def compute_streamline_2D(x, y, z, u, v, w, xy_0): interp = LinearNDInterpolator(xy, uv) def velocity_field_interpolate_reg(xy_, s): - if (xy_ <= -0.5).any() or (xy_ >= 0.5).any(): + if (xy_ < -0.5).any() or (xy_ > 0.5).any(): return np.array([0, 0]) else: return interp(xy_)[0] - #return interp(xy_)[0] - streamline = odeint(velocity_field_interpolate_reg, xy_0, s) # 1000 x 2 + streamline = odeint(velocity_field_interpolate_reg, xy_0, s) # n_points x 2 # Determine point of exit from specified domain if (streamline < -0.5).any(): @@ -431,6 +434,43 @@ def velocity_field_interpolate_reg(xy_, s): return streamline +def shift_points_by_streamline_2D(streamline, x, y, z): + ## I think should be dist by turbine only (not each rotor point); and then + # shift all rotor points by the same amount. + # TODO: How will this work over findices? Need to work that out. + xy = np.concatenate((x.mean(axis=(2,3)),y.mean(axis=(2,3))), axis=0) + + rotor_y = y - y.mean(axis=(2,3), keepdims=True) + rotor_z = z - z.mean(axis=(2,3), keepdims=True) + + # Storage for new points + x_new = x.copy() + y_new = y.copy() + + disps_to_streamline = streamline[:,:,None] - xy[None,:, :] + dists_to_streamline = np.linalg.norm(disps_to_streamline, axis=1) + min_idx = np.argmin(dists_to_streamline, axis=0) + #points = streamline[min_idx, :] + + # At this point, loop over the (downstream) turbines for the shifts + # TODO: for viz/solve for points, are there many more points than turbines? + # If so, maybe need to work on a non-looped version of this. + for tindex in range(len(min_idx)): + streamline_diffs = np.diff(streamline[:min_idx[tindex], :], axis=0) + along_stream_dist = np.sum(np.linalg.norm(streamline_diffs, axis=1)) + off_stream_dist = dists_to_streamline[min_idx[tindex], tindex] + #off_stream_disp = disps_to_streamline[min_idx[tindex], :, tindex] + #disp_x = off_stream_disp[0] + #disp_y = off_stream_disp[1] + + #theta = np.arctan2(disp_y, disp_x) + + x_new[:,tindex,:,:] = streamline[0,0] + along_stream_dist + #y_new[:,tindex,:,:] = streamline[0,1] + disp_x*np.sin(theta) + disp_y*np.cos(theta) + rotor_y[:,tindex,:,:] + y_new[:,tindex,:,:] = streamline[0,1] + off_stream_dist + rotor_y[:,tindex,:,:] + + return x_new, y_new, z + class Loader(yaml.SafeLoader): def __init__(self, stream): From 2fed76fbaa9a7ba89cb25ca077e3474a69835f7c Mon Sep 17 00:00:00 2001 From: misi9170 Date: Wed, 31 Jul 2024 14:49:33 -0600 Subject: [PATCH 18/37] Works over findices, but by looping (may be slow for many findices) --- .../examples_heterogeneous/viz_by_sweep.py | 14 ++-- floris/core/grid.py | 12 +-- floris/utilities.py | 78 ++++++++++++------- 3 files changed, 63 insertions(+), 41 deletions(-) diff --git a/examples/examples_heterogeneous/viz_by_sweep.py b/examples/examples_heterogeneous/viz_by_sweep.py index 5f846871a..9319b777e 100644 --- a/examples/examples_heterogeneous/viz_by_sweep.py +++ b/examples/examples_heterogeneous/viz_by_sweep.py @@ -15,12 +15,12 @@ # # Some wake models may not yet have a visualization method included, for these cases can use # # a slower version which scans a turbine model to produce the horizontal flow -x = np.array([-100, 0, 100, 200, 300, 400, 500, 600]) -v = -np.array([[0, 1, 2, 3, 4, 5, 6, 7]]) +x = np.linspace(-100, 1000, 8) +v = np.linspace(0, 7, 8).reshape(1, 8) het_inflow_config = { "x": np.repeat(x, 2), - "y": np.tile(np.array([-100.0, 100.0]), 8), + "y": np.tile(np.array([-200.0, 200.0]), 8), #"z": np.array([90.0, 90.0, 90.0, 91.0]), "speed_multipliers": 1.0*np.ones((1, 16)), "u": 8.0*np.ones((1, v.shape[1]*2)), @@ -30,8 +30,8 @@ # Set a 2 turbine layout fmodel.set( - layout_x=[0], - layout_y=[0], + layout_x=[0, 500, 0, 500], + layout_y=[0, 0, -100, -100], wind_directions=[270], wind_speeds=[8], turbulence_intensities=[0.06], @@ -42,8 +42,8 @@ fmodel, x_resolution=20, y_resolution=10, - x_bounds=(-100, 500), - y_bounds=(-100, 100), + x_bounds=(-100, 1000), + y_bounds=(-200, 200), ) fig, ax = plt.subplots(figsize=(10, 4)) diff --git a/floris/core/grid.py b/floris/core/grid.py index 6f98add39..6c44844fe 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -262,12 +262,12 @@ def set_grid(self) -> None: self.x_sorted, self.y_sorted, self.z_sorted, - np.array(self.heterogeneous_inflow_config["x"]), - np.array(self.heterogeneous_inflow_config["y"]), - np.array(self.heterogeneous_inflow_config["z"]) if "z" in self.heterogeneous_inflow_config else None, - np.array(self.heterogeneous_inflow_config["u"]), - np.array(self.heterogeneous_inflow_config["v"]), - np.array(self.heterogeneous_inflow_config["w"]) if "w" in self.heterogeneous_inflow_config else None, + self.heterogeneous_inflow_config["x"], + self.heterogeneous_inflow_config["y"], + self.heterogeneous_inflow_config["z"] if "z" in self.heterogeneous_inflow_config else None, + self.heterogeneous_inflow_config["u"], + self.heterogeneous_inflow_config["v"], + self.heterogeneous_inflow_config["w"] if "w" in self.heterogeneous_inflow_config else None, self.n_turbines ) diff --git a/floris/utilities.py b/floris/utilities.py index cb3555b21..045c6e00b 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -302,6 +302,18 @@ def apply_wind_direction_heterogeneity( y_per_turbine = np.repeat(y[:,:,:,:,None], n_turbines, axis=4) z_per_turbine = np.repeat(z[:,:,:,:,None], n_turbines, axis=4) + n_findex = x.shape[0] + + # Ensure all needed elements are arrays + wind_direction_x_points = np.array(wind_direction_x_points) + wind_direction_y_points = np.array(wind_direction_y_points) + wind_direction_u_values = np.array(wind_direction_u_values) + wind_direction_v_values = np.array(wind_direction_v_values) + if wind_direction_z_points is not None: + wind_direction_z_points = np.array(wind_direction_z_points) + if wind_direction_w_values is not None: + wind_direction_w_values = np.array(wind_direction_w_values) + # Compute new relative locations # TEMP________ option = 3 @@ -316,26 +328,33 @@ def apply_wind_direction_heterogeneity( y_per_turbine[0,:,:,:,0] = y_per_turbine[0,:,:,:,0] + 0.2*x_per_turbine[0,:,:,:,0] elif option == 3: # TODO: compute streamlines; convert to point movements. - # TODO: do I need an outer loop over the findices? Maybe; ideally not. - for tindex in range(n_turbines): - # 1. Compute the streamline - streamline = compute_streamline( - wind_direction_x_points, - wind_direction_y_points, - wind_direction_z_points, - wind_direction_u_values, - wind_direction_v_values, - wind_direction_w_values, - np.array([x[0,tindex,:,:].mean(), y[0,tindex,:,:].mean()]) - ) - - # 2. Compute the point movements - x_shifted, y_shifted, z_shifted = shift_points_by_streamline(streamline, x, y, z) - - # 3. Apply to per_turbine values - x_per_turbine[:,:,:,:,tindex] = x_shifted - y_per_turbine[:,:,:,:,tindex] = y_shifted - z_per_turbine[:,:,:,:,tindex] = z_shifted + # TODO: Can I avoid any of these looping operations? Do I need to? + for findex in range(n_findex): + for tindex in range(n_turbines): + # 1. Compute the streamline + streamline = compute_streamline( + wind_direction_x_points, + wind_direction_y_points, + None if wind_direction_z_points is None else wind_direction_z_points[findex,:], + wind_direction_u_values[findex,:], + wind_direction_v_values[findex,:], + None if wind_direction_w_values is None else wind_direction_w_values[findex,:], + np.array([ + x[findex,tindex,:,:].mean(), + y[findex,tindex,:,:].mean(), + z[findex,tindex,:,:].mean() + ]) + ) + + # 2. Compute the point movements + x_shifted, y_shifted, z_shifted = shift_points_by_streamline( + streamline, x[findex,:,:,:], y[findex,:,:,:], z[findex,:,:,:] + ) + + # 3. Apply to per_turbine values + x_per_turbine[findex,:,:,:,tindex] = x_shifted + y_per_turbine[findex,:,:,:,tindex] = y_shifted + z_per_turbine[findex,:,:,:,tindex] = z_shifted else: pass # make no change @@ -380,11 +399,14 @@ def velocity_field_interpolate_reg(xyz_, s): return streamline -def compute_streamline_2D(x, y, z, u, v, w, xy_0, n_points=1000): +def compute_streamline_2D(x, y, z, u, v, w, xyz_0, n_points=1000): """ - Compute streamline starting at xy_0. + Compute streamline starting at xyz_0. z and w will be ignored. """ + + xy_0 = xyz_0[:2] + s = np.linspace(0, 1, n_points) # parametric variable scale = np.array([x.max() - x.min(), y.max() - y.min()]) @@ -393,7 +415,7 @@ def compute_streamline_2D(x, y, z, u, v, w, xy_0, n_points=1000): # Collect xy = np.array([x, y]).T - uv = np.concatenate((u, v), axis=0).T + uv = np.array([u,v]).T # Normalize xy = (xy - offset) / scale @@ -438,10 +460,10 @@ def shift_points_by_streamline_2D(streamline, x, y, z): ## I think should be dist by turbine only (not each rotor point); and then # shift all rotor points by the same amount. # TODO: How will this work over findices? Need to work that out. - xy = np.concatenate((x.mean(axis=(2,3)),y.mean(axis=(2,3))), axis=0) + xy = np.concatenate((x.mean(axis=(1,2))[None,:],y.mean(axis=(1,2))[None,:]), axis=0) - rotor_y = y - y.mean(axis=(2,3), keepdims=True) - rotor_z = z - z.mean(axis=(2,3), keepdims=True) + rotor_y = y - y.mean(axis=(1,2), keepdims=True) + rotor_z = z - z.mean(axis=(1,2), keepdims=True) # Storage for new points x_new = x.copy() @@ -465,9 +487,9 @@ def shift_points_by_streamline_2D(streamline, x, y, z): #theta = np.arctan2(disp_y, disp_x) - x_new[:,tindex,:,:] = streamline[0,0] + along_stream_dist + x_new[tindex,:,:] = streamline[0,0] + along_stream_dist #y_new[:,tindex,:,:] = streamline[0,1] + disp_x*np.sin(theta) + disp_y*np.cos(theta) + rotor_y[:,tindex,:,:] - y_new[:,tindex,:,:] = streamline[0,1] + off_stream_dist + rotor_y[:,tindex,:,:] + y_new[tindex,:,:] = streamline[0,1] + off_stream_dist + rotor_y[tindex,:,:] return x_new, y_new, z From 92b3ae487481ebb15accd299b49ac3715357fd0b Mon Sep 17 00:00:00 2001 From: misi9170 Date: Wed, 31 Jul 2024 17:54:27 -0600 Subject: [PATCH 19/37] Fixes for full flow solve and visualization. --- examples/examples_heterogeneous/bulk_het_4.py | 13 +- examples/examples_heterogeneous/viz_by_viz.py | 58 ++++++++ floris/core/grid.py | 6 + floris/floris_model.py | 10 +- floris/utilities.py | 126 +++++++++--------- 5 files changed, 144 insertions(+), 69 deletions(-) create mode 100644 examples/examples_heterogeneous/viz_by_viz.py diff --git a/examples/examples_heterogeneous/bulk_het_4.py b/examples/examples_heterogeneous/bulk_het_4.py index 326b9e33f..1491496a3 100644 --- a/examples/examples_heterogeneous/bulk_het_4.py +++ b/examples/examples_heterogeneous/bulk_het_4.py @@ -9,7 +9,7 @@ # Get a test fi (single turbine at 0,0) fmodel = FlorisModel("../inputs/gch.yaml") -fmodel.set(layout_x=[0, 0, 600, 600], layout_y=[0, -100, 0, -100]) +fmodel.set(layout_x=[0, 0, 600, 600], layout_y=[0, -300, 0, -300]) fmodel.set( wind_directions=[270.0], @@ -22,8 +22,7 @@ het_inflow_config = { "x": np.array([-1000.0, -1000.0, 1000.0, 1000.0]), - "y": np.array([-1000.0, 1000.0, -1000.0, 1000.0]), - #"z": np.array([90.0, 90.0, 90.0, 91.0]), + "y": np.array([-500.0, 500.0, -500.0, 500.0]), "speed_multipliers": np.array([[1.0, 1.0, 1.0, 1.0]]), "u": np.array([[8.0, 8.0, 8.0, 8.0]]), "v": np.array([[-2.0, 0.0, -2.0, 0.0]]), @@ -40,10 +39,9 @@ fig, ax = plt.subplots(2, 1) fmodel.set( heterogeneous_inflow_config={ - #"wind_directions": [270.0, 270.0, 270.0, 270.0], "x": np.array([-1000.0, -1000.0, 1000.0, 1000.0]), "y": np.array([-1000.0, 1000.0, -1000.0, 1000.0]), - "speed_multipliers":np.array([[1.0, 1.0, 1.0, 1.0]])#np.array([[1.0, 1.2, 1.2, 1.4]]) + "speed_multipliers":np.array([[1.0, 1.0, 1.0, 1.0]]) } ) fmodel.run() @@ -51,6 +49,8 @@ x_resolution=200, y_resolution=100, height=90.0, + x_bounds=(-200, 1000), + y_bounds=(-500, 500), ) visualize_cut_plane( horizontal_plane, @@ -64,7 +64,10 @@ x_resolution=200, y_resolution=100, height=90.0, + x_bounds=(-200, 1000), + y_bounds=(-500, 500), ) + #import ipdb; ipdb.set_trace() visualize_cut_plane( horizontal_plane, ax=ax[1], diff --git a/examples/examples_heterogeneous/viz_by_viz.py b/examples/examples_heterogeneous/viz_by_viz.py new file mode 100644 index 000000000..efa01bbee --- /dev/null +++ b/examples/examples_heterogeneous/viz_by_viz.py @@ -0,0 +1,58 @@ +"""Example: Visualize flow by sweeping turbines + +Demonstrate the use calculate_horizontal_plane_with_turbines + +""" + +import matplotlib.pyplot as plt +import numpy as np + +import floris.flow_visualization as flowviz +from floris import FlorisModel + + +fmodel = FlorisModel("../inputs/gch.yaml") + +# # Some wake models may not yet have a visualization method included, for these cases can use +# # a slower version which scans a turbine model to produce the horizontal flow +x = np.array([-100, 0, 100, 200, 300, 400, 500, 600]) +v = np.array([[0, 1, 2, 3, 4, 5, 6, 7]]) + +het_inflow_config = { + "x": np.repeat(x, 2), + "y": np.tile(np.array([-100.0, 100.0]), 8), + #"z": np.array([90.0, 90.0, 90.0, 91.0]), + "speed_multipliers": 1.0*np.ones((1, 16)), + "u": 8.0*np.ones((1, v.shape[1]*2)), + "v": np.repeat(v, 2, axis=1) + #"v": np.array([[0.0, 0.0, 0.0, 0.0]]), +} + +# Set a 2 turbine layout +fmodel.set( + layout_x=[0, 300], + layout_y=[0, 0], + wind_directions=[270], + wind_speeds=[8], + turbulence_intensities=[0.06], + heterogeneous_inflow_config=het_inflow_config, +) + +horizontal_plane = fmodel.calculate_horizontal_plane( + x_resolution=200, + y_resolution=100, + height=90.0, + x_bounds=(-100, 500), + y_bounds=(-100, 100), +) + +fig, ax = plt.subplots(figsize=(10, 4)) +flowviz.visualize_cut_plane( + horizontal_plane, + ax=ax, + label_contours=True, + title="Horizontal (coarse viz method)", +) + + +plt.show() diff --git a/floris/core/grid.py b/floris/core/grid.py index 6c44844fe..29d192b15 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -259,6 +259,9 @@ def set_grid(self) -> None: # set up new x, y, z coordinates as a mockup self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \ apply_wind_direction_heterogeneity( + x, + y, + z, self.x_sorted, self.y_sorted, self.z_sorted, @@ -636,6 +639,9 @@ def set_grid(self) -> None: # set up new x, y, z coordinates as a mockup self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \ apply_wind_direction_heterogeneity( + x, + y, + z, self.x_sorted, self.y_sorted, self.z_sorted, diff --git a/floris/floris_model.py b/floris/floris_model.py index 3cf8ff8df..469648429 100644 --- a/floris/floris_model.py +++ b/floris/floris_model.py @@ -983,11 +983,13 @@ def set_for_viz(self, findex: int, solver_settings: dict) -> None: self.core.flow_field.heterogeneous_inflow_config[k] ) for k in ["wind_directions", "u", "v", "w"]: # Optional, dependent on findex - # TODO: add findex dimension back in! if k in self.core.flow_field.heterogeneous_inflow_config: - heterogeneous_inflow_config[k] = ( - self.core.flow_field.heterogeneous_inflow_config[k]#[findex:findex+1], - ) + if self.core.flow_field.heterogeneous_inflow_config[k] is not None: + heterogeneous_inflow_config[k] = ( + self.core.flow_field.heterogeneous_inflow_config[k][findex:findex+1,:] + ) + else: + heterogeneous_inflow_config[k] = None else: heterogeneous_inflow_config = None diff --git a/floris/utilities.py b/floris/utilities.py index 045c6e00b..a3a6b5b3c 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -282,9 +282,12 @@ def warp_grid_for_wind_direction_heterogeneity( return x, y_warped, z def apply_wind_direction_heterogeneity( - x, - y, - z, + x_turbines, + y_turbines, + z_turbines, + x_points, + y_points, + z_points, wind_direction_x_points, wind_direction_y_points, wind_direction_z_points, @@ -298,11 +301,11 @@ def apply_wind_direction_heterogeneity( shift_points_by_streamline = shift_points_by_streamline_2D# if wind_direction_z_points is None else shift_points_by_streamline_3D # Initialize storage - x_per_turbine = np.repeat(x[:,:,:,:,None], n_turbines, axis=4) - y_per_turbine = np.repeat(y[:,:,:,:,None], n_turbines, axis=4) - z_per_turbine = np.repeat(z[:,:,:,:,None], n_turbines, axis=4) + x_points_per_turbine = np.repeat(x_points[:,:,:,:,None], n_turbines, axis=4) + y_points_per_turbine = np.repeat(y_points[:,:,:,:,None], n_turbines, axis=4) + z_points_per_turbine = np.repeat(z_points[:,:,:,:,None], n_turbines, axis=4) - n_findex = x.shape[0] + n_findex = x_turbines.shape[0] # Ensure all needed elements are arrays wind_direction_x_points = np.array(wind_direction_x_points) @@ -319,19 +322,23 @@ def apply_wind_direction_heterogeneity( option = 3 if option == 1: # Place turbine 2 directly downstream of turbine 1, as seen by turbine 1 - y_per_turbine[0,2,:,:,1] = y_per_turbine[0,1,:,:,1] + y_points_per_turbine[0,2,:,:,1] = y_points_per_turbine[0,1,:,:,1] # Place turbine 3 out of the direct path of turbine 1, as seen by turbine 1 - y_per_turbine[0,3,:,:,1] = 2*y_per_turbine[0,3,:,:,1] + y_points_per_turbine[0,3,:,:,1] = 2*y_points_per_turbine[0,3,:,:,1] elif option == 2: #y_per_turbine[0,:,:,:,0] = y_per_turbine[0,:,:,:,0] - 0.2*x_per_turbine[0,:,:,:,0] #y_per_turbine[0,:,:,:,1] = y_per_turbine[0,:,:,:,1] + 0.2*x_per_turbine[0,:,:,:,1] - y_per_turbine[0,:,:,:,0] = y_per_turbine[0,:,:,:,0] + 0.2*x_per_turbine[0,:,:,:,0] + y_points_per_turbine[0,:,:,:,0] = y_points_per_turbine[0,:,:,:,0] + 0.2*x_points_per_turbine[0,:,:,:,0] elif option == 3: - # TODO: compute streamlines; convert to point movements. # TODO: Can I avoid any of these looping operations? Do I need to? for findex in range(n_findex): for tindex in range(n_turbines): # 1. Compute the streamline + streamline_start = np.array([ + x_turbines[findex,tindex], + y_turbines[findex,tindex], + z_turbines[findex,tindex] + ]) streamline = compute_streamline( wind_direction_x_points, wind_direction_y_points, @@ -339,65 +346,62 @@ def apply_wind_direction_heterogeneity( wind_direction_u_values[findex,:], wind_direction_v_values[findex,:], None if wind_direction_w_values is None else wind_direction_w_values[findex,:], - np.array([ - x[findex,tindex,:,:].mean(), - y[findex,tindex,:,:].mean(), - z[findex,tindex,:,:].mean() - ]) + streamline_start ) # 2. Compute the point movements x_shifted, y_shifted, z_shifted = shift_points_by_streamline( - streamline, x[findex,:,:,:], y[findex,:,:,:], z[findex,:,:,:] + streamline, x_points[findex,:,:,:], y_points[findex,:,:,:], z_points[findex,:,:,:] ) # 3. Apply to per_turbine values - x_per_turbine[findex,:,:,:,tindex] = x_shifted - y_per_turbine[findex,:,:,:,tindex] = y_shifted - z_per_turbine[findex,:,:,:,tindex] = z_shifted + x_points_per_turbine[findex,:,:,:,tindex] = x_shifted + y_points_per_turbine[findex,:,:,:,tindex] = y_shifted + z_points_per_turbine[findex,:,:,:,tindex] = z_shifted else: pass # make no change # END TEMP____ # Return full set of locations - return x_per_turbine, y_per_turbine, z_per_turbine + return x_points_per_turbine, y_points_per_turbine, z_points_per_turbine -def compute_streamline_3D(x, y, z, u, v, w, xyz_0): - """ - Compute streamline starting at xyz_0. - """ - s = np.linspace(0, 1, 1000) # parametric variable +# def compute_streamline_3D(x, y, z, u, v, w, xyz_0): +# """ +# TODO +# Compute streamline starting at xyz_0. +# """ +# s = np.linspace(0, 1, 1000) # parametric variable - scale = np.array([x.max() - x.min(), y.max() - y.min(), z.max() - z.min()]) - offset = np.array([(x.max() + x.min())/2, (y.max() + y.min())/2, (z.max() + z.min())/2]) - # Need to compute the offset too +# scale = np.array([x.max() - x.min(), y.max() - y.min(), z.max() - z.min()]) +# offset = np.array([(x.max() + x.min())/2, (y.max() + y.min())/2, (z.max() + z.min())/2]) +# # Need to compute the offset too - # Collect - xyz = np.array([x, y, z]).T - uvw = np.concatenate((u, v, w), axis=0).T +# # Collect +# xyz = np.array([x, y, z]).T +# uvw = np.concatenate((u, v, w), axis=0).T - # Normalize - xyz = (xyz - offset) / scale - xyz_0 = (xyz_0 - offset) / scale - # TODO: Loop over findices! Maybe not necessary---can I compute all streamlines at once?? - # Would it help? Currently assuming a single findex. - interp = LinearNDInterpolator(xyz, uvw, fill_value=0.0) +# # Normalize +# xyz = (xyz - offset) / scale +# xyz_0 = (xyz_0 - offset) / scale +# # TODO: Loop over findices! Maybe not necessary---can I compute all streamlines at once?? +# # Would it help? Currently assuming a single findex. +# interp = LinearNDInterpolator(xyz, uvw, fill_value=0.0) - def velocity_field_interpolate_reg(xyz_, s): - # if (xyz_ <= -0.5).any() or (xyz_ >= 0.5).any(): - # print("hmm") - # return np.array([0.0, 0.0, 0.0]) - # else: - # return interp(xyz_)[0] - return interp(xyz_)[0] +# def velocity_field_interpolate_reg(xyz_, s): +# # if (xyz_ <= -0.5).any() or (xyz_ >= 0.5).any(): +# # print("hmm") +# # return np.array([0.0, 0.0, 0.0]) +# # else: +# # return interp(xyz_)[0] +# return interp(xyz_)[0] - streamline = odeint(velocity_field_interpolate_reg, xyz_0, s) # 1000 x 3 +# streamline = odeint(velocity_field_interpolate_reg, xyz_0, s) # 1000 x 3 - # Denormalize - streamline = streamline * scale + offset +# # Denormalize +# streamline = streamline * scale + offset - return streamline +# return streamline def compute_streamline_2D(x, y, z, u, v, w, xyz_0, n_points=1000): """ @@ -462,34 +466,36 @@ def shift_points_by_streamline_2D(streamline, x, y, z): # TODO: How will this work over findices? Need to work that out. xy = np.concatenate((x.mean(axis=(1,2))[None,:],y.mean(axis=(1,2))[None,:]), axis=0) - rotor_y = y - y.mean(axis=(1,2), keepdims=True) + rotor_y = y - y.mean(axis=(1,2), keepdims=True) # TODO: still ok for viz solve? rotor_z = z - z.mean(axis=(1,2), keepdims=True) # Storage for new points x_new = x.copy() y_new = y.copy() - disps_to_streamline = streamline[:,:,None] - xy[None,:, :] + disps_to_streamline = xy[None,:,:] - streamline[:, :,None] dists_to_streamline = np.linalg.norm(disps_to_streamline, axis=1) min_idx = np.argmin(dists_to_streamline, axis=0) #points = streamline[min_idx, :] - # At this point, loop over the (downstream) turbines for the shifts - # TODO: for viz/solve for points, are there many more points than turbines? - # If so, maybe need to work on a non-looped version of this. + # TODO: for viz/solve for points, there are many more points than turbines. + # It would be nice not to have a loop here. for tindex in range(len(min_idx)): streamline_diffs = np.diff(streamline[:min_idx[tindex], :], axis=0) along_stream_dist = np.sum(np.linalg.norm(streamline_diffs, axis=1)) off_stream_dist = dists_to_streamline[min_idx[tindex], tindex] - #off_stream_disp = disps_to_streamline[min_idx[tindex], :, tindex] - #disp_x = off_stream_disp[0] - #disp_y = off_stream_disp[1] + off_stream_disp = disps_to_streamline[min_idx[tindex], :, tindex] + disp_x = off_stream_disp[0] + disp_y = off_stream_disp[1] - #theta = np.arctan2(disp_y, disp_x) + theta = np.arctan2(disp_y, disp_x) x_new[tindex,:,:] = streamline[0,0] + along_stream_dist - #y_new[:,tindex,:,:] = streamline[0,1] + disp_x*np.sin(theta) + disp_y*np.cos(theta) + rotor_y[:,tindex,:,:] - y_new[tindex,:,:] = streamline[0,1] + off_stream_dist + rotor_y[tindex,:,:] + # y_new[tindex,:,:] = streamline[0,1] + disp_x*np.sin(theta) + disp_y*np.cos(theta) + rotor_y[tindex,:,:] + y_new[tindex,:,:] = streamline[0,1] + np.sign(theta)*off_stream_dist + rotor_y[tindex,:,:] + #y_new[tindex,:,:] = streamline[0,1] + off_stream_dist + rotor_y[tindex,:,:] + #if tindex == 1: + # import ipdb; ipdb.set_trace() return x_new, y_new, z From a391e7d949c35005bb50eae79d4794c2e2423d9a Mon Sep 17 00:00:00 2001 From: misi9170 Date: Thu, 1 Aug 2024 09:39:49 -0600 Subject: [PATCH 20/37] Committing first-pass tests. --- tests/wind_direction_heterogeneity_test.py | 128 +++++++++++++++++++++ 1 file changed, 128 insertions(+) create mode 100644 tests/wind_direction_heterogeneity_test.py diff --git a/tests/wind_direction_heterogeneity_test.py b/tests/wind_direction_heterogeneity_test.py new file mode 100644 index 000000000..ba5dd430e --- /dev/null +++ b/tests/wind_direction_heterogeneity_test.py @@ -0,0 +1,128 @@ +from pathlib import Path + +import numpy as np +import pytest + +from floris import FlorisModel +from floris.heterogeneous_map import HeterogeneousMap + + +TEST_DATA = Path(__file__).resolve().parent / "data" +YAML_INPUT = TEST_DATA / "input_full.yaml" + +layout_x = [0.0, 500.0] +layout_y = [0.0, 0.0] + +heterogeneous_wd_inflow_config_2D = { + "x": [-100, -100, 1000, 1000], + "y": [-100, 100, -100, 100], + "u": [[8.0, 8.0, 8.0, 8.0]], + "v": [[0.0, 0.0, 8.0, 8.0]], + "speed_multipliers": [[1.0, 1.0, 1.0, 1.0]] # Currently a necessary input +} + +def test_power(): + fmodel = FlorisModel(configuration=YAML_INPUT) + + # First run without heterogeneity + fmodel.set(layout_x=layout_x, layout_y=layout_y) + fmodel.run() + P_without_het = fmodel.get_turbine_powers()[0,:] + + # Add wind direction heterogeneity and rerun + fmodel.set(heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D) + fmodel.run() + P_with_het = fmodel.get_turbine_powers()[0,:] + + # Confirm upstream power the same; downstream power increased + assert np.allclose(P_without_het[0], P_with_het[0]) + assert P_with_het[1] > P_without_het[1] + +def test_symmetry(): + fmodel = FlorisModel(configuration=YAML_INPUT) + + # Set layout and heterogeneity + fmodel.set( + layout_x=layout_x, + layout_y=layout_y, + heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D + ) + + # Run in this direction + fmodel.run() + P_positive_v = fmodel.get_turbine_powers()[0,:] + + # Switch the sign of the inflow v component and rurun + fmodel.set( + heterogeneous_inflow_config={ + "x": heterogeneous_wd_inflow_config_2D["x"], + "y": heterogeneous_wd_inflow_config_2D["y"], + "u": heterogeneous_wd_inflow_config_2D["u"], + "v": -np.array(heterogeneous_wd_inflow_config_2D["v"]), + "speed_multipliers": heterogeneous_wd_inflow_config_2D["speed_multipliers"] + } + ) + fmodel.run() + P_negative_v = fmodel.get_turbine_powers()[0,:] + + # Confirm symmetry + assert np.allclose(P_positive_v, P_negative_v) + +def test_multiple_findices(): + fmodel = FlorisModel(configuration=YAML_INPUT) + + # Set layout and heterogeneity + fmodel.set( + layout_x=layout_x, + layout_y=layout_y, + heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D + ) + fmodel.run() + P_1findex = fmodel.get_turbine_powers()[0,:] + + # Duplicate the inflow conditons + fmodel.set( + heterogeneous_inflow_config={ + "x": heterogeneous_wd_inflow_config_2D["x"], + "y": heterogeneous_wd_inflow_config_2D["y"], + "u": np.repeat(heterogeneous_wd_inflow_config_2D["u"], 2, axis=0), + "v": np.repeat(heterogeneous_wd_inflow_config_2D["v"], 2, axis=0), + "speed_multipliers": np.repeat( + heterogeneous_wd_inflow_config_2D["speed_multipliers"], 2, axis=0 + ) + }, + wind_directions = [270.0, 270.0], + wind_speeds=[8.0, 8.0], + turbulence_intensities=[0.06, 0.06] + ) + fmodel.run() + P_2findex = fmodel.get_turbine_powers() + + # Confirm the same results + assert np.allclose(P_1findex, P_2findex) + + + + +def test_rotated_wind_direction(): # Fails because turbines are not correctly inside the domain (possibly) + fmodel = FlorisModel(configuration=YAML_INPUT) + fmodel.set( + layout_x=layout_x, + layout_y=layout_y, + heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D + ) + fmodel.run() + P_270deg = fmodel.get_turbine_powers()[0,:] + + # Rotate the wind direction by 90 degrees and rerun + fmodel.set( + wind_directions=[0], + heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D + ) + fmodel.run() + P_0deg = fmodel.get_turbine_powers()[0,:] + + print(P_270deg) + print(P_0deg) + + From d4128534ea58924469d613317cff8b12368d9b68 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Thu, 1 Aug 2024 09:50:15 -0600 Subject: [PATCH 21/37] Remove examples and tests relating to earlier prototype. --- examples/examples_heterogeneous/bulk_het.py | 80 ---------------- examples/examples_heterogeneous/bulk_het_2.py | 49 ---------- examples/examples_heterogeneous/bulk_het_3.py | 71 -------------- tests/bulk_wd_hetergoeneity_unit_test.py | 93 ------------------- 4 files changed, 293 deletions(-) delete mode 100644 examples/examples_heterogeneous/bulk_het.py delete mode 100644 examples/examples_heterogeneous/bulk_het_2.py delete mode 100644 examples/examples_heterogeneous/bulk_het_3.py delete mode 100644 tests/bulk_wd_hetergoeneity_unit_test.py diff --git a/examples/examples_heterogeneous/bulk_het.py b/examples/examples_heterogeneous/bulk_het.py deleted file mode 100644 index 0b5989d4e..000000000 --- a/examples/examples_heterogeneous/bulk_het.py +++ /dev/null @@ -1,80 +0,0 @@ -from pathlib import Path - -import matplotlib.pyplot as plt -import numpy as np - -from floris import FlorisModel - - -# Get a test fi (single turbine at 0,0) -fmodel = FlorisModel("../inputs/gch.yaml") -fmodel.set(layout_x=[0], layout_y=[0]) - -# Directly downstream at 270 degrees -option = False -if not option: - sample_x = [500.0] - sample_y = [0.0] - sample_z = [90.0] - - # Sweep across wind directions - wd_array = np.arange(180, 360, 1) - ws_array = 8.0 * np.ones_like(wd_array) - ti_array = 0.06 * np.ones_like(wd_array) - fmodel.set(wind_directions=wd_array, wind_speeds=ws_array, turbulence_intensities=ti_array) - - # Standard case; expect minimum to be at 270 degrees - u_at_points_0 = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z) - - # Now, apply bulk wind direction heterogeneity to the flow - fmodel.set( - heterogeneous_inflow_config={ - "bulk_wd_change": [[0.0, 10.0]]*180, # -10 degree change, CW - "bulk_wd_x": [[0, 500.0]]*180 - } - ) # TODO: Build something that checks the dimensions of the inputs here - u_at_points_1 = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z) - - fmodel.set( - heterogeneous_inflow_config={ - "bulk_wd_change": [[0.0, -10.0]]*180, # -10 degree change, CW - "bulk_wd_x": [[0, 500.0]]*180 - } - ) # TODO: Build something that checks the dimensions of the inputs here - u_at_points_2 = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z) - - - fig, ax = plt.subplots(1,1) - ax.plot(wd_array, u_at_points_0, label="Standard", color="black") - ax.plot(wd_array, u_at_points_1, label="Shifted 10") - ax.plot(wd_array, u_at_points_2, label="Shifted -10") - -else: - ### Let's try something else - sample_x = [500.0]*100 - sample_y = np.linspace(-500, 500, 100) - sample_z = [90.0]*100 - - fmodel.set( - wind_directions=[270.0], - wind_speeds=[8.0], - turbulence_intensities=[0.06], - heterogeneous_inflow_config={} - ) - u_at_points_00 = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z).flatten() - - fmodel.set( - heterogeneous_inflow_config={ - "bulk_wd_change": [[0.0, -30.0]], # -10 degree change, CW - "bulk_wd_x": [[0, 500.0]] - }, - ) - u_at_points_01 = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z).flatten() - - fig, ax = plt.subplots(1,1) - ax.plot(sample_y, u_at_points_00, label="Standard", color="black") - ax.plot(sample_y, u_at_points_01, label="Shifted 10") - -plt.show() - - diff --git a/examples/examples_heterogeneous/bulk_het_2.py b/examples/examples_heterogeneous/bulk_het_2.py deleted file mode 100644 index da1f49f1d..000000000 --- a/examples/examples_heterogeneous/bulk_het_2.py +++ /dev/null @@ -1,49 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np - -import floris.layout_visualization as layoutviz -from floris import FlorisModel -from floris.flow_visualization import visualize_cut_plane - - -# Get a test fi (single turbine at 0,0) -fmodel = FlorisModel("../inputs/gch.yaml") -fmodel.set(layout_x=[0, 600], layout_y=[0, 0]) - -# Directly downstream at 270 degrees -### Let's try something else -sample_x = [500.0]*100 -sample_y = np.linspace(-500, 500, 100) -sample_z = [90.0]*100 - -fmodel.set( - wind_directions=[270.0], - wind_speeds=[8.0], - turbulence_intensities=[0.06], - heterogeneous_inflow_config={ - "bulk_wd_change": [[0.0, -20.0, 20.0, -20.0]], # -10 degree change, CW - #"bulk_wd_change": [[0.0, 0.0, 0.0, 0.0]], # -10 degree change, CCW - "bulk_wd_x": [[0, 600, 1200, 1800]] - }, -) -horizontal_plane = fmodel.calculate_horizontal_plane( - x_resolution=200, - y_resolution=100, - height=90.0, -) - -# Plot the flow field with rotors -fig, ax = plt.subplots() -visualize_cut_plane( - horizontal_plane, - ax=ax, - label_contours=False, - #title="Horizontal Flow with Turbine Rotors and labels", -) - -# Plot the turbine rotors -ax.grid() - -plt.show() - - diff --git a/examples/examples_heterogeneous/bulk_het_3.py b/examples/examples_heterogeneous/bulk_het_3.py deleted file mode 100644 index c9e00b80e..000000000 --- a/examples/examples_heterogeneous/bulk_het_3.py +++ /dev/null @@ -1,71 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np - -import floris.layout_visualization as layoutviz -from floris import FlorisModel -from floris.flow_visualization import visualize_cut_plane - -visualize = True - -# Get a test fi (single turbine at 0,0) -fmodel = FlorisModel("../inputs/gch.yaml") -fmodel.set(layout_x=[0, 0, 600, 600], layout_y=[0, -100, 0, -100]) - -fmodel.set( - wind_directions=[270.0], - wind_speeds=[8.0], - turbulence_intensities=[0.06], -) -fmodel.run() -P_wo_het = fmodel.get_turbine_powers()/1e6 - -het_inflow_config = { - "wind_directions": [270.0, 280.0, 280.0, 290.0], - "x": np.array([-1000.0, -1000.0, 1000.0, 1000.0]), - "y": np.array([-1000.0, 1000.0, -1000.0, 1000.0]), - "speed_multipliers":np.array([[1.0, 1.0, 1.0, 1.0]])#np.array([[1.0, 1.2, 1.2, 1.4]]) -} - -fmodel.set(heterogeneous_inflow_config=het_inflow_config) -fmodel.run() -P_w_het = fmodel.get_turbine_powers()/1e6 - -print("Difference (MW):", P_w_het - P_wo_het) - -if visualize: - fig, ax = plt.subplots(2, 1) - fmodel.set( - heterogeneous_inflow_config={ - #"wind_directions": [270.0, 270.0, 270.0, 270.0], - "x": np.array([-1000.0, -1000.0, 1000.0, 1000.0]), - "y": np.array([-1000.0, 1000.0, -1000.0, 1000.0]), - "speed_multipliers":np.array([[1.0, 1.0, 1.0, 1.0]])#np.array([[1.0, 1.2, 1.2, 1.4]]) - } - ) - fmodel.run() - horizontal_plane = fmodel.calculate_horizontal_plane( - x_resolution=200, - y_resolution=100, - height=90.0, - ) - visualize_cut_plane( - horizontal_plane, - ax=ax[0], - label_contours=False, - title="Without WD het" - ) - - fmodel.set(heterogeneous_inflow_config=het_inflow_config) - horizontal_plane = fmodel.calculate_horizontal_plane( - x_resolution=200, - y_resolution=100, - height=90.0, - ) - visualize_cut_plane( - horizontal_plane, - ax=ax[1], - label_contours=False, - title="With WD het" - ) - - plt.show() diff --git a/tests/bulk_wd_hetergoeneity_unit_test.py b/tests/bulk_wd_hetergoeneity_unit_test.py deleted file mode 100644 index 3d8721062..000000000 --- a/tests/bulk_wd_hetergoeneity_unit_test.py +++ /dev/null @@ -1,93 +0,0 @@ -import logging -from pathlib import Path - -import numpy as np -import pytest - -from floris import ( - FlorisModel, - TimeSeries, - WindRose, -) -from floris.optimization.layout_optimization.layout_optimization_base import ( - LayoutOptimization, -) -from floris.optimization.layout_optimization.layout_optimization_scipy import ( - LayoutOptimizationScipy, -) -from floris.wind_data import WindDataBase - - -TEST_DATA = Path(__file__).resolve().parent / "data" -YAML_INPUT = TEST_DATA / "input_full.yaml" - -def test_bulk_wd_heterogeneity_turbine(): - fmodel = FlorisModel(configuration=YAML_INPUT) - fmodel.set(layout_x=[0, 600], layout_y=[0, 0]) - - # Run straight as a benchmark - fmodel.set(wind_directions=[270.0], wind_speeds=[8.0], turbulence_intensities=[0.06]) - fmodel.run() - powers_straight = fmodel.get_turbine_powers()[0,:] - - - # Single wind direction, with two wind direction shifts as well as straight - fmodel.set( - wind_directions=[270.0, 270.0, 270.0], - wind_speeds=[8.0, 8.0, 8.0], - turbulence_intensities=[0.06, 0.06, 0.06], - heterogeneous_inflow_config={ - "bulk_wd_change": [[0.0, 0.0], [0.0, -10.0], [0.0, 10.0]], # 10 degree changes - "bulk_wd_x": [[0, 1000], [0, 1000], [0, 1000]] # Past downstream turbine - } - ) - - fmodel.run() - powers = fmodel.get_turbine_powers() - - assert (powers_straight == powers[0,:]).all() # Verify straight case - - assert (powers[:,0] == powers[0,0]).all() # Upstream turbine not affected - assert powers[0,1] < powers[1,1] # wake shifted away from downstream turbine - assert np.allclose(powers[1,1], powers[2,1]) # Shifted wake is symmetric - - # Check works for other directions too - fmodel.set( - wind_directions=[225.0, 225.0, 225.0], - layout_x=[0, 600/np.sqrt(2)], - layout_y=[0, 600/np.sqrt(2)], - ) - fmodel.run() - powers_2 = fmodel.get_turbine_powers() - - assert np.allclose(powers_2, powers) - -def test_bulk_wd_heterogeneity_flow_field(): - # Get a test fi (single turbine at 0,0) - fmodel = FlorisModel(configuration=YAML_INPUT) - - # Directly downstream at 270 degrees - sample_x = [500.0] - sample_y = [0.0] - sample_z = [90.0] - - # Sweep across wind directions - wd_array = np.arange(180, 360, 1) - ws_array = 8.0 * np.ones_like(wd_array) - ti_array = 0.06 * np.ones_like(wd_array) - fmodel.set(wind_directions=wd_array, wind_speeds=ws_array, turbulence_intensities=ti_array) - - # Standard case; expect minimum to be at 270 degrees - u_at_points = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z) - assert (wd_array[np.argmin(u_at_points)] == 270) - - # Now, apply bulk wind direction heterogeneity to the flow - fmodel.set( - heterogeneous_inflow_config={ - "bulk_wd_change": [[0.0, -10.0]]*180, # -10 degree change, CW - "bulk_wd_x": [[0, 500.0]]*180 - } - ) # TODO: Build something that checks the dimensions of the inputs here - - u_at_points = fmodel.sample_flow_at_points(sample_x, sample_y, sample_z) - assert (wd_array[np.argmin(u_at_points)] == 280) From e34ec8521c80dc25b2146c365809d5bb57c58e7b Mon Sep 17 00:00:00 2001 From: misi9170 Date: Thu, 1 Aug 2024 10:43:40 -0600 Subject: [PATCH 22/37] Remove running code for earlier prototype. --- floris/core/flow_field.py | 7 ----- floris/core/grid.py | 40 +++------------------------ floris/floris_model.py | 2 +- floris/utilities.py | 57 +-------------------------------------- 4 files changed, 5 insertions(+), 101 deletions(-) diff --git a/floris/core/flow_field.py b/floris/core/flow_field.py index f9e3bab2a..ce2d6529b 100644 --- a/floris/core/flow_field.py +++ b/floris/core/flow_field.py @@ -312,13 +312,6 @@ def _speed_heterogeneity(self): else: return False - @property - def _bulk_wd_heterogeneity(self): - if self.heterogeneous_inflow_config is not None: - return "bulk_wd_change" in self.heterogeneous_inflow_config - else: - return False - @staticmethod def interpolate_multiplier_xy(x: NDArrayFloat, y: NDArrayFloat, diff --git a/floris/core/grid.py b/floris/core/grid.py index cfeba192f..6f3f337fa 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -18,8 +18,7 @@ from floris.utilities import ( reverse_rotate_coordinates_rel_west, rotate_coordinates_rel_west, - warp_grid_for_wind_direction_heterogeneity, - apply_wind_direction_heterogeneity, + apply_wind_direction_heterogeneity_warping, ) @@ -258,7 +257,7 @@ def set_grid(self) -> None: self.use_turbine_specific_layouts = True # set up new x, y, z coordinates as a mockup self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \ - apply_wind_direction_heterogeneity( + apply_wind_direction_heterogeneity_warping( x, y, z, @@ -321,13 +320,6 @@ def set_grid(self) -> None: self.wind_directions, self.turbine_coordinates ) - if (self.heterogeneous_inflow_config is not None - and "bulk_wd_change" in self.heterogeneous_inflow_config): - wd_het_x_points = np.array(self.heterogeneous_inflow_config["bulk_wd_x"]) - wd_het_values = np.array(self.heterogeneous_inflow_config["bulk_wd_change"]) - x, y, z = warp_grid_for_wind_direction_heterogeneity( - x, y, z, wd_het_x_points + x.min(), wd_het_values - ) # Coefficients cubature_coefficients = TurbineCubatureGrid.get_cubature_coefficients(self.grid_resolution) @@ -510,13 +502,6 @@ def set_grid(self) -> None: self.wind_directions, self.turbine_coordinates ) - if (self.heterogeneous_inflow_config is not None - and "bulk_wd_change" in self.heterogeneous_inflow_config): - wd_het_x_points = np.array(self.heterogeneous_inflow_config["bulk_wd_x"]) - wd_het_values = np.array(self.heterogeneous_inflow_config["bulk_wd_change"]) - x, y, z = warp_grid_for_wind_direction_heterogeneity( - x, y, z, wd_het_x_points + x.min(), wd_het_values - ) # Construct the arrays storing the grid points eps = 0.01 @@ -649,7 +634,7 @@ def set_grid(self) -> None: self.use_turbine_specific_layouts = True # set up new x, y, z coordinates as a mockup self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \ - apply_wind_direction_heterogeneity( + apply_wind_direction_heterogeneity_warping( x, y, z, @@ -675,18 +660,6 @@ def set_grid(self) -> None: y_center_of_rotation=self.y_center_of_rotation, ) - if (self.heterogeneous_inflow_config is not None - and "bulk_wd_change" in self.heterogeneous_inflow_config): - wd_het_x_points = np.array(self.heterogeneous_inflow_config["bulk_wd_x"]) - wd_het_values = np.array(self.heterogeneous_inflow_config["bulk_wd_change"]) - x_points, y_points, z_points = warp_grid_for_wind_direction_heterogeneity( - x_points, y_points, z_points, wd_het_x_points + x.min(), wd_het_values - ) - - self.x_sorted = x_points[None, :, :, :] - self.y_sorted = y_points[None, :, :, :] - self.z_sorted = z_points[None, :, :, :] - @define class PointsGrid(Grid): """ @@ -733,13 +706,6 @@ def set_grid(self) -> None: x_center_of_rotation=self.x_center_of_rotation, y_center_of_rotation=self.y_center_of_rotation ) - if (self.heterogeneous_inflow_config is not None - and "bulk_wd_change" in self.heterogeneous_inflow_config): - wd_het_x_points = np.array(self.heterogeneous_inflow_config["bulk_wd_x"]) - wd_het_values = np.array(self.heterogeneous_inflow_config["bulk_wd_change"]) - x, y, z = warp_grid_for_wind_direction_heterogeneity( - x, y, z, wd_het_x_points, wd_het_values - ) self.x_sorted = x[:,:,None,None] self.y_sorted = y[:,:,None,None] self.z_sorted = z[:,:,None,None] diff --git a/floris/floris_model.py b/floris/floris_model.py index 70d4c9ce4..233ee0ba6 100644 --- a/floris/floris_model.py +++ b/floris/floris_model.py @@ -1049,7 +1049,7 @@ def set_for_viz(self, findex: int, solver_settings: dict) -> None: heterogeneous_inflow_config[k] = ( self.core.flow_field.heterogeneous_inflow_config[k] ) - for k in ["wind_directions", "u", "v", "w"]: # Optional, dependent on findex + for k in ["u", "v", "w"]: # Optional, dependent on findex if k in self.core.flow_field.heterogeneous_inflow_config: if self.core.flow_field.heterogeneous_inflow_config[k] is not None: heterogeneous_inflow_config[k] = ( diff --git a/floris/utilities.py b/floris/utilities.py index 5b2417db8..0e764522f 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -327,62 +327,7 @@ def reverse_rotate_coordinates_rel_west( return grid_x_reversed, grid_y_reversed, grid_z_reversed -def warp_grid_for_wind_direction_heterogeneity( - x, - y, - z, - wd_het_x_points, - wd_het_values, - ): - - y_warped = y.copy() - - for i in range(wd_het_values.shape[1]-1): - - x0 = wd_het_x_points[:,i:i+1] - x1 = wd_het_x_points[:,i+1:i+2] - - # Compute radius of curvature for the warp (not sign swap in tand) - warp_radius = (x1 - x0) / tand(wd_het_values[:,i:i+1] - wd_het_values[:,i+1:i+2]) - - ### Below is code for a "true" warp in x and y, but is not complete. - # # Compute the angle of warping for each grid point (radians) - # warp_theta = np.arctan2(x - x0, np.abs(warp_radius))*np.sign(warp_radius) - - # # Apply warp - # import ipdb; ipdb.set_trace() - # x_warped = np.where( # Also add end points - # (x >= x0) & (x <= x1) & (np.abs(warp_radius) < np.inf), - # warp_radius * np.sin(warp_theta) + y * np.sin(warp_theta) + x0, - # x - # ) - # y_warped = np.where( - # (x >= x0) & (x <= x1) & (np.abs(warp_radius) < np.inf), - # warp_radius * np.cos(warp_theta) + y * np.cos(warp_theta) - warp_radius, - # y - # ) - - ### Instead, use a simpler y "slide" approach - if len(x.shape) == len(x1.shape): - x1_rep = np.repeat(x1, x.shape[1], axis=1) # Usual solve - else: - x1_rep = x1[:,:,None] # Full field solve - x_mod = np.where( - x > x1_rep, - x1_rep, - x - ) - y_warped = np.where( - (np.abs(warp_radius) != np.inf ) & (x >= x0), - (-warp_radius + np.sign(warp_radius) * np.sqrt(warp_radius**2 - (x_mod - x0)**2) - + y_warped), - y_warped - ) - - # Return warped grid - return x, y_warped, z - -def apply_wind_direction_heterogeneity( +def apply_wind_direction_heterogeneity_warping( x_turbines, y_turbines, z_turbines, From bb38c5ce5b54e0c6d867ef8332754c1703444e3e Mon Sep 17 00:00:00 2001 From: misi9170 Date: Thu, 1 Aug 2024 11:04:48 -0600 Subject: [PATCH 23/37] Consolidate and rename examples. --- ...py => 00x_wind_direction_heterogeneity.py} | 0 ..._direction_heterogeneity_visualization.py} | 30 +++++++--- .../examples_heterogeneous/viz_by_sweep.py | 58 ------------------- 3 files changed, 22 insertions(+), 66 deletions(-) rename examples/examples_heterogeneous/{bulk_het_4.py => 00x_wind_direction_heterogeneity.py} (100%) rename examples/examples_heterogeneous/{viz_by_viz.py => 00x_wind_direction_heterogeneity_visualization.py} (63%) delete mode 100644 examples/examples_heterogeneous/viz_by_sweep.py diff --git a/examples/examples_heterogeneous/bulk_het_4.py b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity.py similarity index 100% rename from examples/examples_heterogeneous/bulk_het_4.py rename to examples/examples_heterogeneous/00x_wind_direction_heterogeneity.py diff --git a/examples/examples_heterogeneous/viz_by_viz.py b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py similarity index 63% rename from examples/examples_heterogeneous/viz_by_viz.py rename to examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py index efa01bbee..a1402234c 100644 --- a/examples/examples_heterogeneous/viz_by_viz.py +++ b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py @@ -13,6 +13,9 @@ fmodel = FlorisModel("../inputs/gch.yaml") +viz_by_sweep = False # TODO: sweep not working correctly. +# Maybe be to do with rotor points, or how calculate_horizontal_plane_with_turbines is implemented + # # Some wake models may not yet have a visualization method included, for these cases can use # # a slower version which scans a turbine model to produce the horizontal flow x = np.array([-100, 0, 100, 200, 300, 400, 500, 600]) @@ -38,20 +41,31 @@ heterogeneous_inflow_config=het_inflow_config, ) -horizontal_plane = fmodel.calculate_horizontal_plane( - x_resolution=200, - y_resolution=100, - height=90.0, - x_bounds=(-100, 500), - y_bounds=(-100, 100), -) +if viz_by_sweep: + horizontal_plane = flowviz.calculate_horizontal_plane_with_turbines( + fmodel, + x_resolution=20, + y_resolution=10, + x_bounds=(-100, 500), + y_bounds=(-100, 100), + ) + title = "Coarse turbine scan method" +else: + horizontal_plane = fmodel.calculate_horizontal_plane( + x_resolution=200, + y_resolution=100, + height=90.0, + x_bounds=(-100, 500), + y_bounds=(-100, 100), + ) + title = "Standard flow visualization calculation" fig, ax = plt.subplots(figsize=(10, 4)) flowviz.visualize_cut_plane( horizontal_plane, ax=ax, label_contours=True, - title="Horizontal (coarse viz method)", + title=title ) diff --git a/examples/examples_heterogeneous/viz_by_sweep.py b/examples/examples_heterogeneous/viz_by_sweep.py deleted file mode 100644 index 9319b777e..000000000 --- a/examples/examples_heterogeneous/viz_by_sweep.py +++ /dev/null @@ -1,58 +0,0 @@ -"""Example: Visualize flow by sweeping turbines - -Demonstrate the use calculate_horizontal_plane_with_turbines - -""" - -import matplotlib.pyplot as plt -import numpy as np - -import floris.flow_visualization as flowviz -from floris import FlorisModel - - -fmodel = FlorisModel("../inputs/gch.yaml") - -# # Some wake models may not yet have a visualization method included, for these cases can use -# # a slower version which scans a turbine model to produce the horizontal flow -x = np.linspace(-100, 1000, 8) -v = np.linspace(0, 7, 8).reshape(1, 8) - -het_inflow_config = { - "x": np.repeat(x, 2), - "y": np.tile(np.array([-200.0, 200.0]), 8), - #"z": np.array([90.0, 90.0, 90.0, 91.0]), - "speed_multipliers": 1.0*np.ones((1, 16)), - "u": 8.0*np.ones((1, v.shape[1]*2)), - "v": np.repeat(v, 2, axis=1) - #"v": np.array([[0.0, 0.0, 0.0, 0.0]]), -} - -# Set a 2 turbine layout -fmodel.set( - layout_x=[0, 500, 0, 500], - layout_y=[0, 0, -100, -100], - wind_directions=[270], - wind_speeds=[8], - turbulence_intensities=[0.06], - heterogeneous_inflow_config=het_inflow_config, -) - -horizontal_plane_scan_turbine = flowviz.calculate_horizontal_plane_with_turbines( - fmodel, - x_resolution=20, - y_resolution=10, - x_bounds=(-100, 1000), - y_bounds=(-200, 200), -) - -fig, ax = plt.subplots(figsize=(10, 4)) -flowviz.visualize_cut_plane( - horizontal_plane_scan_turbine, - ax=ax, - label_contours=True, - title="Horizontal (coarse turbine scan method)", -) - - -plt.show() From a428113c65b42587c0d2c914722d1afc4489f7ca Mon Sep 17 00:00:00 2001 From: misi9170 Date: Thu, 1 Aug 2024 14:18:17 -0600 Subject: [PATCH 24/37] Change turbine ordering so that viz by sweep works again. --- ...d_direction_heterogeneity_visualization.py | 3 +-- floris/core/grid.py | 17 ++++++------ floris/utilities.py | 26 ++++++++++++++++--- 3 files changed, 32 insertions(+), 14 deletions(-) diff --git a/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py index a1402234c..54663169b 100644 --- a/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py +++ b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py @@ -13,8 +13,7 @@ fmodel = FlorisModel("../inputs/gch.yaml") -viz_by_sweep = False # TODO: sweep not working correctly. -# Maybe be to do with rotor points, or how calculate_horizontal_plane_with_turbines is implemented +viz_by_sweep = False # # Some wake models may not yet have a visualization method included, for these cases can use # # a slower version which scans a turbine model to produce the horizontal flow diff --git a/floris/core/grid.py b/floris/core/grid.py index 6f3f337fa..e7b5731a6 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -258,9 +258,9 @@ def set_grid(self) -> None: # set up new x, y, z coordinates as a mockup self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \ apply_wind_direction_heterogeneity_warping( - x, - y, - z, + self.x_sorted.mean(axis=(2,3)), + self.y_sorted.mean(axis=(2,3)), + self.z_sorted.mean(axis=(2,3)), self.x_sorted, self.y_sorted, self.z_sorted, @@ -270,7 +270,6 @@ def set_grid(self) -> None: self.heterogeneous_inflow_config["u"], self.heterogeneous_inflow_config["v"], self.heterogeneous_inflow_config["w"] if "w" in self.heterogeneous_inflow_config else None, - self.n_turbines ) # Now calculate grid coordinates in original frame (from 270 deg perspective) @@ -632,12 +631,13 @@ def set_grid(self) -> None: if (self.heterogeneous_inflow_config is not None and "u" in self.heterogeneous_inflow_config): self.use_turbine_specific_layouts = True - # set up new x, y, z coordinates as a mockup + + # Apply wind direction heterogeneity to generate turbine-specific layouts self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \ apply_wind_direction_heterogeneity_warping( - x, - y, - z, + np.take_along_axis(x, x.argsort(axis=1), axis=1), + np.take_along_axis(y, x.argsort(axis=1), axis=1), + np.take_along_axis(z, x.argsort(axis=1), axis=1), self.x_sorted, self.y_sorted, self.z_sorted, @@ -647,7 +647,6 @@ def set_grid(self) -> None: self.heterogeneous_inflow_config["u"], self.heterogeneous_inflow_config["v"], self.heterogeneous_inflow_config["w"] if "w" in self.heterogeneous_inflow_config else None, - self.n_turbines ) # Now calculate grid coordinates in original frame (from 270 deg perspective) self.x_sorted_inertial_frame, self.y_sorted_inertial_frame, self.z_sorted_inertial_frame = \ diff --git a/floris/utilities.py b/floris/utilities.py index 0e764522f..276cbe87d 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -340,19 +340,34 @@ def apply_wind_direction_heterogeneity_warping( wind_direction_u_values, wind_direction_v_values, wind_direction_w_values, - n_turbines ): + """ + Args: + x_turbines (np.ndarray): x-coordinates of the turbines (n_findex x n_turbines). + y_turbines (np.ndarray): y-coordinates of the turbines (n_findex x n_turbines). + z_turbines (np.ndarray): z-coordinates of the turbines (n_findex x n_turbines). + x_points (np.ndarray): x-coordinates to be warped (n_findex x n_turbines x n_grid x n_grid). + y_points (np.ndarray): y-coordinates to be warped (n_findex x n_turbines x n_grid x n_grid). + z_points (np.ndarray): z-coordinates to be warped (n_findex x n_turbines x n_grid x n_grid). + wind_direction_x_points (np.ndarray): x-coordinates of the specified flow points (n_findex x n_points). + wind_direction_y_points (np.ndarray): y-coordinates of the specified flow points (n_findex x n_points). + wind_direction_z_points (np.ndarray): z-coordinates of the specified flow points (n_findex x n_points). + wind_direction_u_values (np.ndarray): u-velocity values at the specified flow points (n_findex x n_points). + wind_direction_v_values (np.ndarray): v-velocity values at the specified flow points (n_findex x n_points). + wind_direction_w_values (np.ndarray): w-velocity values at the specified flow points (n_findex x n_points). + n_turbines (int): Number of turbines. + """ compute_streamline = compute_streamline_2D# if wind_direction_z_points is None else compute_streamline_3D shift_points_by_streamline = shift_points_by_streamline_2D# if wind_direction_z_points is None else shift_points_by_streamline_3D + n_findex, n_turbines = x_turbines.shape + # Initialize storage x_points_per_turbine = np.repeat(x_points[:,:,:,:,None], n_turbines, axis=4) y_points_per_turbine = np.repeat(y_points[:,:,:,:,None], n_turbines, axis=4) z_points_per_turbine = np.repeat(z_points[:,:,:,:,None], n_turbines, axis=4) - n_findex = x_turbines.shape[0] - # Ensure all needed elements are arrays wind_direction_x_points = np.array(wind_direction_x_points) wind_direction_y_points = np.array(wind_direction_y_points) @@ -385,6 +400,11 @@ def apply_wind_direction_heterogeneity_warping( y_turbines[findex,tindex], z_turbines[findex,tindex] ]) + # streamline_start = np.array([ + # x_points[findex,tindex,:,:].mean(), + # y_points[findex,tindex,:,:].mean(), + # z_points[findex,tindex,:,:].mean(), + # ]) streamline = compute_streamline( wind_direction_x_points, wind_direction_y_points, From 64916e4bdf0d3671331e49ddf134bd928439f5ad Mon Sep 17 00:00:00 2001 From: misi9170 Date: Thu, 1 Aug 2024 14:20:10 -0600 Subject: [PATCH 25/37] Remove unneeded debugging code. --- floris/utilities.py | 77 +++++++++++++++++---------------------------- 1 file changed, 28 insertions(+), 49 deletions(-) diff --git a/floris/utilities.py b/floris/utilities.py index 276cbe87d..78aef8bc6 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -379,55 +379,34 @@ def apply_wind_direction_heterogeneity_warping( wind_direction_w_values = np.array(wind_direction_w_values) # Compute new relative locations - # TEMP________ - option = 3 - if option == 1: - # Place turbine 2 directly downstream of turbine 1, as seen by turbine 1 - y_points_per_turbine[0,2,:,:,1] = y_points_per_turbine[0,1,:,:,1] - # Place turbine 3 out of the direct path of turbine 1, as seen by turbine 1 - y_points_per_turbine[0,3,:,:,1] = 2*y_points_per_turbine[0,3,:,:,1] - elif option == 2: - #y_per_turbine[0,:,:,:,0] = y_per_turbine[0,:,:,:,0] - 0.2*x_per_turbine[0,:,:,:,0] - #y_per_turbine[0,:,:,:,1] = y_per_turbine[0,:,:,:,1] + 0.2*x_per_turbine[0,:,:,:,1] - y_points_per_turbine[0,:,:,:,0] = y_points_per_turbine[0,:,:,:,0] + 0.2*x_points_per_turbine[0,:,:,:,0] - elif option == 3: - # TODO: Can I avoid any of these looping operations? Do I need to? - for findex in range(n_findex): - for tindex in range(n_turbines): - # 1. Compute the streamline - streamline_start = np.array([ - x_turbines[findex,tindex], - y_turbines[findex,tindex], - z_turbines[findex,tindex] - ]) - # streamline_start = np.array([ - # x_points[findex,tindex,:,:].mean(), - # y_points[findex,tindex,:,:].mean(), - # z_points[findex,tindex,:,:].mean(), - # ]) - streamline = compute_streamline( - wind_direction_x_points, - wind_direction_y_points, - None if wind_direction_z_points is None else wind_direction_z_points[findex,:], - wind_direction_u_values[findex,:], - wind_direction_v_values[findex,:], - None if wind_direction_w_values is None else wind_direction_w_values[findex,:], - streamline_start - ) - - # 2. Compute the point movements - x_shifted, y_shifted, z_shifted = shift_points_by_streamline( - streamline, x_points[findex,:,:,:], y_points[findex,:,:,:], z_points[findex,:,:,:] - ) - - # 3. Apply to per_turbine values - x_points_per_turbine[findex,:,:,:,tindex] = x_shifted - y_points_per_turbine[findex,:,:,:,tindex] = y_shifted - z_points_per_turbine[findex,:,:,:,tindex] = z_shifted - else: - pass # make no change - - # END TEMP____ + # TODO: Can I avoid any of these looping operations? Do I need to? + for findex in range(n_findex): + for tindex in range(n_turbines): + # 1. Compute the streamline + streamline_start = np.array([ + x_turbines[findex,tindex], + y_turbines[findex,tindex], + z_turbines[findex,tindex] + ]) + streamline = compute_streamline( + wind_direction_x_points, + wind_direction_y_points, + None if wind_direction_z_points is None else wind_direction_z_points[findex,:], + wind_direction_u_values[findex,:], + wind_direction_v_values[findex,:], + None if wind_direction_w_values is None else wind_direction_w_values[findex,:], + streamline_start + ) + + # 2. Compute the point movements + x_shifted, y_shifted, z_shifted = shift_points_by_streamline( + streamline, x_points[findex,:,:,:], y_points[findex,:,:,:], z_points[findex,:,:,:] + ) + + # 3. Apply to per_turbine values + x_points_per_turbine[findex,:,:,:,tindex] = x_shifted + y_points_per_turbine[findex,:,:,:,tindex] = y_shifted + z_points_per_turbine[findex,:,:,:,tindex] = z_shifted # Return full set of locations return x_points_per_turbine, y_points_per_turbine, z_points_per_turbine From 049b7defd5eaca5586314345608f90c1baae07b5 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Thu, 1 Aug 2024 14:53:38 -0600 Subject: [PATCH 26/37] Require that Grid children be instantiated with required keyword args so that Grid() can have heterogeneous_inflow_config take a default of None. --- floris/core/grid.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/floris/core/grid.py b/floris/core/grid.py index e7b5731a6..27c41cbec 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -53,7 +53,7 @@ class Grid(ABC, BaseClass): turbine_diameters: NDArrayFloat = field(converter=floris_array_converter) wind_directions: NDArrayFloat = field(converter=floris_array_converter) grid_resolution: int | Iterable = field() - heterogeneous_inflow_config: dict | None = field() + heterogeneous_inflow_config: dict | None = field(init=True, default=None) n_turbines: int = field(init=False) n_findex: int = field(init=False) @@ -111,7 +111,7 @@ def grid_resolution_validator(self, instance: attrs.Attribute, value: int | Iter def set_grid(self) -> None: raise NotImplementedError("Grid.set_grid") -@define +@define(kw_only=True) class TurbineGrid(Grid): """See `Grid` for more details. @@ -283,7 +283,7 @@ def set_grid(self) -> None: y_center_of_rotation=self.y_center_of_rotation, ) -@define +@define(kw_only=True) class TurbineCubatureGrid(Grid): """ This grid type arranges points throughout the swept area of the rotor based on the cubature @@ -464,7 +464,7 @@ def get_cubature_coefficients(cls, N: int): "B": np.pi/N, } -@define +@define(kw_only=True) class FlowFieldGrid(Grid): """ Args: @@ -533,7 +533,7 @@ def set_grid(self) -> None: y_center_of_rotation=self.y_center_of_rotation, ) -@define +@define(kw_only=True) class FlowFieldPlanarGrid(Grid): """ Args: @@ -659,7 +659,7 @@ def set_grid(self) -> None: y_center_of_rotation=self.y_center_of_rotation, ) -@define +@define(kw_only=True) class PointsGrid(Grid): """ Args: From 7d1c8ed24d217a446969d72eccc238e4822861ae Mon Sep 17 00:00:00 2001 From: misi9170 Date: Fri, 2 Aug 2024 13:16:32 -0600 Subject: [PATCH 27/37] Formatting. --- .../00x_wind_direction_heterogeneity.py | 3 +- ...d_direction_heterogeneity_visualization.py | 2 +- floris/core/grid.py | 16 ++++--- floris/utilities.py | 48 +++++++++++-------- tests/wind_direction_heterogeneity_test.py | 6 +-- 5 files changed, 44 insertions(+), 31 deletions(-) diff --git a/examples/examples_heterogeneous/00x_wind_direction_heterogeneity.py b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity.py index 1491496a3..9fe69d335 100644 --- a/examples/examples_heterogeneous/00x_wind_direction_heterogeneity.py +++ b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity.py @@ -5,6 +5,7 @@ from floris import FlorisModel from floris.flow_visualization import visualize_cut_plane + visualize = True # Get a test fi (single turbine at 0,0) @@ -36,7 +37,7 @@ print("Difference (MW):", P_w_het - P_wo_het) if visualize: - fig, ax = plt.subplots(2, 1) + fig, ax = plt.subplots(2, 1, figsize=(4,8)) fmodel.set( heterogeneous_inflow_config={ "x": np.array([-1000.0, -1000.0, 1000.0, 1000.0]), diff --git a/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py index 54663169b..fb0d61b6d 100644 --- a/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py +++ b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py @@ -13,7 +13,7 @@ fmodel = FlorisModel("../inputs/gch.yaml") -viz_by_sweep = False +viz_by_sweep = True # # Some wake models may not yet have a visualization method included, for these cases can use # # a slower version which scans a turbine model to produce the horizontal flow diff --git a/floris/core/grid.py b/floris/core/grid.py index 27c41cbec..aac9f1b66 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -16,9 +16,9 @@ NDArrayInt, ) from floris.utilities import ( + apply_wind_direction_heterogeneity_warping, reverse_rotate_coordinates_rel_west, rotate_coordinates_rel_west, - apply_wind_direction_heterogeneity_warping, ) @@ -266,10 +266,12 @@ def set_grid(self) -> None: self.z_sorted, self.heterogeneous_inflow_config["x"], self.heterogeneous_inflow_config["y"], - self.heterogeneous_inflow_config["z"] if "z" in self.heterogeneous_inflow_config else None, + self.heterogeneous_inflow_config["z"] + if "z" in self.heterogeneous_inflow_config else None, self.heterogeneous_inflow_config["u"], self.heterogeneous_inflow_config["v"], - self.heterogeneous_inflow_config["w"] if "w" in self.heterogeneous_inflow_config else None, + self.heterogeneous_inflow_config["w"] + if "w" in self.heterogeneous_inflow_config else None, ) # Now calculate grid coordinates in original frame (from 270 deg perspective) @@ -631,7 +633,7 @@ def set_grid(self) -> None: if (self.heterogeneous_inflow_config is not None and "u" in self.heterogeneous_inflow_config): self.use_turbine_specific_layouts = True - + # Apply wind direction heterogeneity to generate turbine-specific layouts self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \ apply_wind_direction_heterogeneity_warping( @@ -643,10 +645,12 @@ def set_grid(self) -> None: self.z_sorted, self.heterogeneous_inflow_config["x"], self.heterogeneous_inflow_config["y"], - self.heterogeneous_inflow_config["z"] if "z" in self.heterogeneous_inflow_config else None, + self.heterogeneous_inflow_config["z"] + if "z" in self.heterogeneous_inflow_config else None, self.heterogeneous_inflow_config["u"], self.heterogeneous_inflow_config["v"], - self.heterogeneous_inflow_config["w"] if "w" in self.heterogeneous_inflow_config else None, + self.heterogeneous_inflow_config["w"] + if "w" in self.heterogeneous_inflow_config else None, ) # Now calculate grid coordinates in original frame (from 270 deg perspective) self.x_sorted_inertial_frame, self.y_sorted_inertial_frame, self.z_sorted_inertial_frame = \ diff --git a/floris/utilities.py b/floris/utilities.py index 78aef8bc6..4ab90057c 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -14,8 +14,8 @@ import numpy as np import yaml from attrs import define, field -from scipy.interpolate import LinearNDInterpolator from scipy.integrate import odeint +from scipy.interpolate import LinearNDInterpolator from floris.type_dec import floris_array_converter, NDArrayFloat @@ -349,17 +349,24 @@ def apply_wind_direction_heterogeneity_warping( x_points (np.ndarray): x-coordinates to be warped (n_findex x n_turbines x n_grid x n_grid). y_points (np.ndarray): y-coordinates to be warped (n_findex x n_turbines x n_grid x n_grid). z_points (np.ndarray): z-coordinates to be warped (n_findex x n_turbines x n_grid x n_grid). - wind_direction_x_points (np.ndarray): x-coordinates of the specified flow points (n_findex x n_points). - wind_direction_y_points (np.ndarray): y-coordinates of the specified flow points (n_findex x n_points). - wind_direction_z_points (np.ndarray): z-coordinates of the specified flow points (n_findex x n_points). - wind_direction_u_values (np.ndarray): u-velocity values at the specified flow points (n_findex x n_points). - wind_direction_v_values (np.ndarray): v-velocity values at the specified flow points (n_findex x n_points). - wind_direction_w_values (np.ndarray): w-velocity values at the specified flow points (n_findex x n_points). - n_turbines (int): Number of turbines. + wind_direction_x_points (np.ndarray): x-coordinates of the specified flow points + (n_findex x n_points). + wind_direction_y_points (np.ndarray): y-coordinates of the specified flow points + (n_findex x n_points). + wind_direction_z_points (np.ndarray): z-coordinates of the specified flow points + (n_findex x n_points). + wind_direction_u_values (np.ndarray): u-velocity values at the specified flow points + (n_findex x n_points). + wind_direction_v_values (np.ndarray): v-velocity values at the specified flow points + (n_findex x n_points). + wind_direction_w_values (np.ndarray): w-velocity values at the specified flow points + (n_findex x n_points). """ - compute_streamline = compute_streamline_2D# if wind_direction_z_points is None else compute_streamline_3D - shift_points_by_streamline = shift_points_by_streamline_2D# if wind_direction_z_points is None else shift_points_by_streamline_3D + # TODO: create switch that uses 3D versions, possible, based on whether + # wind_direction_z_points is None + compute_streamline = compute_streamline_2D + shift_points_by_streamline = shift_points_by_streamline_2D n_findex, n_turbines = x_turbines.shape @@ -402,7 +409,7 @@ def apply_wind_direction_heterogeneity_warping( x_shifted, y_shifted, z_shifted = shift_points_by_streamline( streamline, x_points[findex,:,:,:], y_points[findex,:,:,:], z_points[findex,:,:,:] ) - + # 3. Apply to per_turbine values x_points_per_turbine[findex,:,:,:,tindex] = x_shifted y_points_per_turbine[findex,:,:,:,tindex] = y_shifted @@ -496,7 +503,7 @@ def velocity_field_interpolate_reg(xy_, s): raise ValueError("Streamline could not be propagated inside the domain.") if exit_index == 1000: raise ValueError("Streamline did not reach the edge of the domain.") - + # Copy points after exit streamline[exit_index:] = streamline[exit_index-1] @@ -506,13 +513,13 @@ def velocity_field_interpolate_reg(xy_, s): return streamline def shift_points_by_streamline_2D(streamline, x, y, z): - ## I think should be dist by turbine only (not each rotor point); and then + ## I think should be dist by turbine only (not each rotor point); and then # shift all rotor points by the same amount. # TODO: How will this work over findices? Need to work that out. xy = np.concatenate((x.mean(axis=(1,2))[None,:],y.mean(axis=(1,2))[None,:]), axis=0) rotor_y = y - y.mean(axis=(1,2), keepdims=True) # TODO: still ok for viz solve? - rotor_z = z - z.mean(axis=(1,2), keepdims=True) + z - z.mean(axis=(1,2), keepdims=True) # Storage for new points x_new = x.copy() @@ -532,11 +539,12 @@ def shift_points_by_streamline_2D(streamline, x, y, z): off_stream_disp = disps_to_streamline[min_idx[tindex], :, tindex] disp_x = off_stream_disp[0] disp_y = off_stream_disp[1] - + theta = np.arctan2(disp_y, disp_x) - + x_new[tindex,:,:] = streamline[0,0] + along_stream_dist - # y_new[tindex,:,:] = streamline[0,1] + disp_x*np.sin(theta) + disp_y*np.cos(theta) + rotor_y[tindex,:,:] + # y_new[tindex,:,:] = streamline[0,1] + disp_x*np.sin(theta) + disp_y*np.cos(theta)\ + # + rotor_y[tindex,:,:] y_new[tindex,:,:] = streamline[0,1] + np.sign(theta)*off_stream_dist + rotor_y[tindex,:,:] #y_new[tindex,:,:] = streamline[0,1] + off_stream_dist + rotor_y[tindex,:,:] #if tindex == 1: @@ -668,14 +676,14 @@ def update_model_args( Update the x_sorted, y_sorted, and z_sorted entries of the model_args dict. """ # Check that new x, y, z are the same length as the old - if (x.shape != model_args["x"].shape or + if (x.shape != model_args["x"].shape or y.shape != model_args["y"].shape or z.shape != model_args["z"].shape): raise ValueError("Size of new x,y,z arrays must match existing ones.") - + # Update the x, y, z values model_args["x"] = x model_args["y"] = y model_args["z"] = z - + return model_args diff --git a/tests/wind_direction_heterogeneity_test.py b/tests/wind_direction_heterogeneity_test.py index ba5dd430e..144ced765 100644 --- a/tests/wind_direction_heterogeneity_test.py +++ b/tests/wind_direction_heterogeneity_test.py @@ -37,7 +37,7 @@ def test_power(): # Confirm upstream power the same; downstream power increased assert np.allclose(P_without_het[0], P_with_het[0]) assert P_with_het[1] > P_without_het[1] - + def test_symmetry(): fmodel = FlorisModel(configuration=YAML_INPUT) @@ -103,8 +103,8 @@ def test_multiple_findices(): - -def test_rotated_wind_direction(): # Fails because turbines are not correctly inside the domain (possibly) +# I think failing because turbines are not correctly inside the domain +def test_rotated_wind_direction(): fmodel = FlorisModel(configuration=YAML_INPUT) fmodel.set( layout_x=layout_x, From 1fc4c32faa0b2aa97e28a62c35f5d2f3b6e898dd Mon Sep 17 00:00:00 2001 From: misi9170 Date: Mon, 5 Aug 2024 12:37:24 -0600 Subject: [PATCH 28/37] Working on both 3d and rotations --- floris/core/grid.py | 48 ++++++++++++++++++++++++++++++++++++++++++--- floris/utilities.py | 43 +++++++++++++++------------------------- 2 files changed, 61 insertions(+), 30 deletions(-) diff --git a/floris/core/grid.py b/floris/core/grid.py index aac9f1b66..e84798424 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -255,7 +255,51 @@ def set_grid(self) -> None: if (self.heterogeneous_inflow_config is not None and "u" in self.heterogeneous_inflow_config): self.use_turbine_specific_layouts = True - # set up new x, y, z coordinates as a mockup + + # Expand z dimension if not provided + # TODO: Is there already a validator somewhere else that this could be done? + if self.heterogeneous_inflow_config["z"] is None: + n_het_points = len(self.heterogeneous_inflow_config["x"]) + self.heterogeneous_inflow_config["z"] = np.concatenate( + (np.zeros(n_het_points), 1000*np.ones(n_het_points)) + ) # Assumes no turbines above 1000m + + # Repeat all other elements + self.heterogeneous_inflow_config["x"] = np.tile( + self.heterogeneous_inflow_config["x"], + 2 + ) + self.heterogeneous_inflow_config["y"] = np.tile( + self.heterogeneous_inflow_config["y"], + 2 + ) + self.heterogeneous_inflow_config["u"] = np.tile( + self.heterogeneous_inflow_config["u"], + (1, 2) + ) + self.heterogeneous_inflow_config["v"] = np.tile( + self.heterogeneous_inflow_config["v"], + (1, 2) + ) + + + # Rotate the input x, y, z to match the wind direction + het_coords = np.concatenate( + ( + np.array(self.heterogeneous_inflow_config["x"])[:, None], + np.array(self.heterogeneous_inflow_config["y"])[:, None], + np.array(self.heterogeneous_inflow_config["z"])[:, None] + ), + axis=1 + ) + x_het, y_het, z_het, _, _ = rotate_coordinates_rel_west( + self.wind_directions, + het_coords, + x_center_of_rotation=self.x_center_of_rotation, + y_center_of_rotation=self.y_center_of_rotation + ) + + self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \ apply_wind_direction_heterogeneity_warping( self.x_sorted.mean(axis=(2,3)), @@ -270,8 +314,6 @@ def set_grid(self) -> None: if "z" in self.heterogeneous_inflow_config else None, self.heterogeneous_inflow_config["u"], self.heterogeneous_inflow_config["v"], - self.heterogeneous_inflow_config["w"] - if "w" in self.heterogeneous_inflow_config else None, ) # Now calculate grid coordinates in original frame (from 270 deg perspective) diff --git a/floris/utilities.py b/floris/utilities.py index 4ab90057c..7b507c662 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -339,7 +339,6 @@ def apply_wind_direction_heterogeneity_warping( wind_direction_z_points, wind_direction_u_values, wind_direction_v_values, - wind_direction_w_values, ): """ Args: @@ -365,7 +364,7 @@ def apply_wind_direction_heterogeneity_warping( # TODO: create switch that uses 3D versions, possible, based on whether # wind_direction_z_points is None - compute_streamline = compute_streamline_2D + #compute_streamline = compute_streamline_2D shift_points_by_streamline = shift_points_by_streamline_2D n_findex, n_turbines = x_turbines.shape @@ -376,14 +375,8 @@ def apply_wind_direction_heterogeneity_warping( z_points_per_turbine = np.repeat(z_points[:,:,:,:,None], n_turbines, axis=4) # Ensure all needed elements are arrays - wind_direction_x_points = np.array(wind_direction_x_points) - wind_direction_y_points = np.array(wind_direction_y_points) wind_direction_u_values = np.array(wind_direction_u_values) wind_direction_v_values = np.array(wind_direction_v_values) - if wind_direction_z_points is not None: - wind_direction_z_points = np.array(wind_direction_z_points) - if wind_direction_w_values is not None: - wind_direction_w_values = np.array(wind_direction_w_values) # Compute new relative locations # TODO: Can I avoid any of these looping operations? Do I need to? @@ -398,10 +391,9 @@ def apply_wind_direction_heterogeneity_warping( streamline = compute_streamline( wind_direction_x_points, wind_direction_y_points, - None if wind_direction_z_points is None else wind_direction_z_points[findex,:], + wind_direction_z_points, wind_direction_u_values[findex,:], wind_direction_v_values[findex,:], - None if wind_direction_w_values is None else wind_direction_w_values[findex,:], streamline_start ) @@ -455,39 +447,36 @@ def apply_wind_direction_heterogeneity_warping( # return streamline -def compute_streamline_2D(x, y, z, u, v, w, xyz_0, n_points=1000): +def compute_streamline(x, y, z, u, v, xyz_0, n_points=1000): """ Compute streamline starting at xyz_0. - z and w will be ignored. """ - xy_0 = xyz_0[:2] - s = np.linspace(0, 1, n_points) # parametric variable - scale = np.array([x.max() - x.min(), y.max() - y.min()]) - offset = np.array([(x.max() + x.min())/2, (y.max() + y.min())/2]) + scale = np.array([x.max() - x.min(), y.max() - y.min(), z.max() - z.min()]) + offset = np.array([(x.max() + x.min())/2, (y.max() + y.min())/2, (z.max() + z.min())/2]) # Need to compute the offset too # Collect - xy = np.array([x, y]).T - uv = np.array([u,v]).T + xyz = np.array([x, y, z]).T + uv = np.array([u, v, np.zeros_like(u)]).T # streamlines evolve in (x,y) space # Normalize - xy = (xy - offset) / scale - xy_0 = (xy_0 - offset) / scale + xyz = (xyz - offset) / scale + xyz_0 = (xyz_0 - offset) / scale # TODO: Loop over findices! Maybe not necessary---can I compute all streamlines at once?? # Would it help? Currently assuming a single findex. - interp = LinearNDInterpolator(xy, uv) + interp = LinearNDInterpolator(xyz, uv) - def velocity_field_interpolate_reg(xy_, s): - if (xy_ < -0.5).any() or (xy_ > 0.5).any(): - return np.array([0, 0]) + def velocity_field_interpolate_reg(xyz_, s): + if (xyz_ < -0.5).any() or (xyz_ > 0.5).any(): + return np.array([0, 0, 0]) else: - return interp(xy_)[0] + return interp(xyz_)[0] - streamline = odeint(velocity_field_interpolate_reg, xy_0, s) # n_points x 2 + streamline = odeint(velocity_field_interpolate_reg, xyz_0, s) # n_points x 3 # Determine point of exit from specified domain if (streamline < -0.5).any(): @@ -512,7 +501,7 @@ def velocity_field_interpolate_reg(xy_, s): return streamline -def shift_points_by_streamline_2D(streamline, x, y, z): +def shift_points_by_streamline(streamline, x, y, z): ## I think should be dist by turbine only (not each rotor point); and then # shift all rotor points by the same amount. # TODO: How will this work over findices? Need to work that out. From d4823449a8b4148ce2def5b4ca26d3f8d18c805a Mon Sep 17 00:00:00 2001 From: misi9170 Date: Mon, 5 Aug 2024 12:50:33 -0600 Subject: [PATCH 29/37] Remove z capabilities; working on rotation. --- floris/core/grid.py | 62 ++++++++++++--------------------- floris/utilities.py | 84 ++++++++++----------------------------------- 2 files changed, 41 insertions(+), 105 deletions(-) diff --git a/floris/core/grid.py b/floris/core/grid.py index e84798424..95a84211a 100644 --- a/floris/core/grid.py +++ b/floris/core/grid.py @@ -256,50 +256,22 @@ def set_grid(self) -> None: and "u" in self.heterogeneous_inflow_config): self.use_turbine_specific_layouts = True - # Expand z dimension if not provided - # TODO: Is there already a validator somewhere else that this could be done? - if self.heterogeneous_inflow_config["z"] is None: - n_het_points = len(self.heterogeneous_inflow_config["x"]) - self.heterogeneous_inflow_config["z"] = np.concatenate( - (np.zeros(n_het_points), 1000*np.ones(n_het_points)) - ) # Assumes no turbines above 1000m - - # Repeat all other elements - self.heterogeneous_inflow_config["x"] = np.tile( - self.heterogeneous_inflow_config["x"], - 2 - ) - self.heterogeneous_inflow_config["y"] = np.tile( - self.heterogeneous_inflow_config["y"], - 2 - ) - self.heterogeneous_inflow_config["u"] = np.tile( - self.heterogeneous_inflow_config["u"], - (1, 2) - ) - self.heterogeneous_inflow_config["v"] = np.tile( - self.heterogeneous_inflow_config["v"], - (1, 2) - ) - - # Rotate the input x, y, z to match the wind direction het_coords = np.concatenate( ( np.array(self.heterogeneous_inflow_config["x"])[:, None], np.array(self.heterogeneous_inflow_config["y"])[:, None], - np.array(self.heterogeneous_inflow_config["z"])[:, None] + np.zeros((len(self.heterogeneous_inflow_config["x"]), 1)) ), axis=1 ) - x_het, y_het, z_het, _, _ = rotate_coordinates_rel_west( + x_het, y_het, _, _, _ = rotate_coordinates_rel_west( self.wind_directions, het_coords, x_center_of_rotation=self.x_center_of_rotation, y_center_of_rotation=self.y_center_of_rotation ) - self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \ apply_wind_direction_heterogeneity_warping( self.x_sorted.mean(axis=(2,3)), @@ -308,10 +280,8 @@ def set_grid(self) -> None: self.x_sorted, self.y_sorted, self.z_sorted, - self.heterogeneous_inflow_config["x"], - self.heterogeneous_inflow_config["y"], - self.heterogeneous_inflow_config["z"] - if "z" in self.heterogeneous_inflow_config else None, + x_het, + y_het, self.heterogeneous_inflow_config["u"], self.heterogeneous_inflow_config["v"], ) @@ -676,6 +646,22 @@ def set_grid(self) -> None: and "u" in self.heterogeneous_inflow_config): self.use_turbine_specific_layouts = True + # Rotate the input x, y, z to match the wind direction + het_coords = np.concatenate( + ( + np.array(self.heterogeneous_inflow_config["x"])[:, None], + np.array(self.heterogeneous_inflow_config["y"])[:, None], + np.zeros((len(self.heterogeneous_inflow_config["x"]), 1)) + ), + axis=1 + ) + x_het, y_het, _, _, _ = rotate_coordinates_rel_west( + self.wind_directions, + het_coords, + x_center_of_rotation=self.x_center_of_rotation, + y_center_of_rotation=self.y_center_of_rotation + ) + # Apply wind direction heterogeneity to generate turbine-specific layouts self.x_sorted_per_turbine, self.y_sorted_per_turbine, self.z_sorted_per_turbine = \ apply_wind_direction_heterogeneity_warping( @@ -685,14 +671,10 @@ def set_grid(self) -> None: self.x_sorted, self.y_sorted, self.z_sorted, - self.heterogeneous_inflow_config["x"], - self.heterogeneous_inflow_config["y"], - self.heterogeneous_inflow_config["z"] - if "z" in self.heterogeneous_inflow_config else None, + x_het, + y_het, self.heterogeneous_inflow_config["u"], self.heterogeneous_inflow_config["v"], - self.heterogeneous_inflow_config["w"] - if "w" in self.heterogeneous_inflow_config else None, ) # Now calculate grid coordinates in original frame (from 270 deg perspective) self.x_sorted_inertial_frame, self.y_sorted_inertial_frame, self.z_sorted_inertial_frame = \ diff --git a/floris/utilities.py b/floris/utilities.py index 7b507c662..bfc5b43a6 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -336,7 +336,6 @@ def apply_wind_direction_heterogeneity_warping( z_points, wind_direction_x_points, wind_direction_y_points, - wind_direction_z_points, wind_direction_u_values, wind_direction_v_values, ): @@ -352,8 +351,6 @@ def apply_wind_direction_heterogeneity_warping( (n_findex x n_points). wind_direction_y_points (np.ndarray): y-coordinates of the specified flow points (n_findex x n_points). - wind_direction_z_points (np.ndarray): z-coordinates of the specified flow points - (n_findex x n_points). wind_direction_u_values (np.ndarray): u-velocity values at the specified flow points (n_findex x n_points). wind_direction_v_values (np.ndarray): v-velocity values at the specified flow points @@ -364,9 +361,6 @@ def apply_wind_direction_heterogeneity_warping( # TODO: create switch that uses 3D versions, possible, based on whether # wind_direction_z_points is None - #compute_streamline = compute_streamline_2D - shift_points_by_streamline = shift_points_by_streamline_2D - n_findex, n_turbines = x_turbines.shape # Initialize storage @@ -385,13 +379,11 @@ def apply_wind_direction_heterogeneity_warping( # 1. Compute the streamline streamline_start = np.array([ x_turbines[findex,tindex], - y_turbines[findex,tindex], - z_turbines[findex,tindex] + y_turbines[findex,tindex] ]) - streamline = compute_streamline( - wind_direction_x_points, - wind_direction_y_points, - wind_direction_z_points, + streamline = compute_2D_streamline( + wind_direction_x_points[findex,:], + wind_direction_y_points[findex,:], wind_direction_u_values[findex,:], wind_direction_v_values[findex,:], streamline_start @@ -410,73 +402,36 @@ def apply_wind_direction_heterogeneity_warping( # Return full set of locations return x_points_per_turbine, y_points_per_turbine, z_points_per_turbine -# def compute_streamline_3D(x, y, z, u, v, w, xyz_0): -# """ -# TODO -# Compute streamline starting at xyz_0. -# """ -# s = np.linspace(0, 1, 1000) # parametric variable - -# scale = np.array([x.max() - x.min(), y.max() - y.min(), z.max() - z.min()]) -# offset = np.array([(x.max() + x.min())/2, (y.max() + y.min())/2, (z.max() + z.min())/2]) -# # Need to compute the offset too - -# # Collect -# xyz = np.array([x, y, z]).T -# uvw = np.concatenate((u, v, w), axis=0).T - -# # Normalize -# xyz = (xyz - offset) / scale -# xyz_0 = (xyz_0 - offset) / scale -# # TODO: Loop over findices! Maybe not necessary---can I compute all streamlines at once?? -# # Would it help? Currently assuming a single findex. -# interp = LinearNDInterpolator(xyz, uvw, fill_value=0.0) - -# def velocity_field_interpolate_reg(xyz_, s): -# # if (xyz_ <= -0.5).any() or (xyz_ >= 0.5).any(): -# # print("hmm") -# # return np.array([0.0, 0.0, 0.0]) -# # else: -# # return interp(xyz_)[0] -# return interp(xyz_)[0] - -# streamline = odeint(velocity_field_interpolate_reg, xyz_0, s) # 1000 x 3 - -# # Denormalize -# streamline = streamline * scale + offset - -# return streamline - -def compute_streamline(x, y, z, u, v, xyz_0, n_points=1000): +def compute_2D_streamline(x, y, u, v, xy_0, n_points=1000): """ - Compute streamline starting at xyz_0. + Compute 2D streamline (in x-y plane) starting at xy_0, according to (u,v) + velocity field. """ - s = np.linspace(0, 1, n_points) # parametric variable - scale = np.array([x.max() - x.min(), y.max() - y.min(), z.max() - z.min()]) - offset = np.array([(x.max() + x.min())/2, (y.max() + y.min())/2, (z.max() + z.min())/2]) + scale = np.array([x.max() - x.min(), y.max() - y.min()]) + offset = np.array([(x.max() + x.min())/2, (y.max() + y.min())/2]) # Need to compute the offset too # Collect - xyz = np.array([x, y, z]).T - uv = np.array([u, v, np.zeros_like(u)]).T # streamlines evolve in (x,y) space + xy = np.array([x, y]).T + uv = np.array([u,v]).T # Normalize - xyz = (xyz - offset) / scale - xyz_0 = (xyz_0 - offset) / scale + xy = (xy - offset) / scale + xy_0 = (xy_0 - offset) / scale # TODO: Loop over findices! Maybe not necessary---can I compute all streamlines at once?? # Would it help? Currently assuming a single findex. - interp = LinearNDInterpolator(xyz, uv) + interp = LinearNDInterpolator(xy, uv) - def velocity_field_interpolate_reg(xyz_, s): - if (xyz_ < -0.5).any() or (xyz_ > 0.5).any(): - return np.array([0, 0, 0]) + def velocity_field_interpolate_reg(xy_, s): + if (xy_ < -0.5).any() or (xy_ > 0.5).any(): + return np.array([0, 0]) else: - return interp(xyz_)[0] + return interp(xy_)[0] - streamline = odeint(velocity_field_interpolate_reg, xyz_0, s) # n_points x 3 + streamline = odeint(velocity_field_interpolate_reg, xy_0, s) # n_points x 2 # Determine point of exit from specified domain if (streamline < -0.5).any(): @@ -508,7 +463,6 @@ def shift_points_by_streamline(streamline, x, y, z): xy = np.concatenate((x.mean(axis=(1,2))[None,:],y.mean(axis=(1,2))[None,:]), axis=0) rotor_y = y - y.mean(axis=(1,2), keepdims=True) # TODO: still ok for viz solve? - z - z.mean(axis=(1,2), keepdims=True) # Storage for new points x_new = x.copy() From 796e14685a56b382d931acdb75314300d1ef0c42 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Mon, 5 Aug 2024 14:31:51 -0600 Subject: [PATCH 30/37] Set up rotation tests correctly; now passing. --- tests/wind_direction_heterogeneity_test.py | 121 +++++++++++++++++---- 1 file changed, 100 insertions(+), 21 deletions(-) diff --git a/tests/wind_direction_heterogeneity_test.py b/tests/wind_direction_heterogeneity_test.py index 144ced765..366493893 100644 --- a/tests/wind_direction_heterogeneity_test.py +++ b/tests/wind_direction_heterogeneity_test.py @@ -10,14 +10,14 @@ TEST_DATA = Path(__file__).resolve().parent / "data" YAML_INPUT = TEST_DATA / "input_full.yaml" -layout_x = [0.0, 500.0] -layout_y = [0.0, 0.0] +layout_x = np.array([0.0, 500.0]) +layout_y = np.array([0.0, 0.0]) heterogeneous_wd_inflow_config_2D = { - "x": [-100, -100, 1000, 1000], - "y": [-100, 100, -100, 100], + "x": [-1000, -1000, 1000, 1000], + "y": [-1000, 1000, -1000, 1000], "u": [[8.0, 8.0, 8.0, 8.0]], - "v": [[0.0, 0.0, 8.0, 8.0]], + "v": [[0.0, 0.0, 4.0, 4.0]], "speed_multipliers": [[1.0, 1.0, 1.0, 1.0]] # Currently a necessary input } @@ -68,6 +68,30 @@ def test_symmetry(): # Confirm symmetry assert np.allclose(P_positive_v, P_negative_v) +def test_input_types(): + fmodel = FlorisModel(configuration=YAML_INPUT) + + # First run without heterogeneity + fmodel.set( + layout_x=layout_x, + layout_y=layout_y, + heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D + ) + fmodel.run() + P_lists = fmodel.get_turbine_powers() + + # Convert all elements of dictionary to arrays + for k in heterogeneous_wd_inflow_config_2D: + heterogeneous_wd_inflow_config_2D[k] = np.array(heterogeneous_wd_inflow_config_2D[k]) + assert isinstance(heterogeneous_wd_inflow_config_2D["u"], np.ndarray) + + fmodel.set(heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D) + fmodel.run() + P_arrays = fmodel.get_turbine_powers() + + # Confirm upstream power the same + assert np.allclose(P_lists, P_arrays) + def test_multiple_findices(): fmodel = FlorisModel(configuration=YAML_INPUT) @@ -101,28 +125,83 @@ def test_multiple_findices(): # Confirm the same results assert np.allclose(P_1findex, P_2findex) - - -# I think failing because turbines are not correctly inside the domain -def test_rotated_wind_direction(): +def test_flipped_direction(): fmodel = FlorisModel(configuration=YAML_INPUT) + + # Set up a case with 2 findices, central turbines, aligned in fmodel.set( - layout_x=layout_x, - layout_y=layout_y, - heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D + wind_directions=[270, 90], + wind_speeds=[8.0, 8.0], + turbulence_intensities=[0.06, 0.06], + layout_x=np.array([-250.0, 250.0]), + layout_y=np.array([0.0, 0.0]), + heterogeneous_inflow_config={ + "x": heterogeneous_wd_inflow_config_2D["x"], + "y": heterogeneous_wd_inflow_config_2D["y"], + "u": np.concatenate( + ( + heterogeneous_wd_inflow_config_2D["u"], + np.flip(heterogeneous_wd_inflow_config_2D["u"], axis=1) + ), + axis=0 + ), + "v": np.concatenate( + ( + heterogeneous_wd_inflow_config_2D["v"], + np.flip(heterogeneous_wd_inflow_config_2D["v"], axis=1) + ), + axis=0 + ), + "speed_multipliers": np.repeat( + heterogeneous_wd_inflow_config_2D["speed_multipliers"], 2, axis=0 + ) + } ) fmodel.run() - P_270deg = fmodel.get_turbine_powers()[0,:] + powers = fmodel.get_turbine_powers() + assert np.allclose(powers[0,:], np.flip(powers[:,1])) + +def test_rotational_symmetry(): + layout_x_symmetry = np.array([-250.0, 250.0]) + layout_y_symmetry = np.array([0.0, 0.0]) - # Rotate the wind direction by 90 degrees and rerun + fmodel = FlorisModel(configuration=YAML_INPUT) fmodel.set( - wind_directions=[0], + layout_x=layout_x_symmetry, + layout_y=layout_y_symmetry, heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D ) fmodel.run() - P_0deg = fmodel.get_turbine_powers()[0,:] - - print(P_270deg) - print(P_0deg) - - + P_base = fmodel.get_turbine_powers() + + # Rotate to a series of new positions, and rotate layout also + wind_directions_test = np.linspace(0, 360, 10) + + for wd in wind_directions_test: + het_wd_inflow_test = heterogeneous_wd_inflow_config_2D.copy() + het_wd_inflow_test["x"] = \ + np.cos(np.deg2rad(270-wd))*np.array(heterogeneous_wd_inflow_config_2D["x"]) \ + - np.sin(np.deg2rad(270-wd))*np.array(heterogeneous_wd_inflow_config_2D["y"]) + het_wd_inflow_test["y"] = \ + np.sin(np.deg2rad(270-wd))*np.array(heterogeneous_wd_inflow_config_2D["x"]) \ + + np.cos(np.deg2rad(270-wd))*np.array(heterogeneous_wd_inflow_config_2D["y"]) + + layout_x_test = ( + np.cos(np.deg2rad(270-wd))*layout_x_symmetry + - np.sin(np.deg2rad(270-wd))*layout_y_symmetry + ) + layout_y_test = ( + np.sin(np.deg2rad(270-wd))*layout_x_symmetry + + np.cos(np.deg2rad(270-wd))*layout_y_symmetry + ) + + fmodel.set( + wind_directions=[wd], + layout_x=layout_x_test, + layout_y=layout_y_test, + heterogeneous_inflow_config=het_wd_inflow_test + ) + fmodel.run() + P_new = fmodel.get_turbine_powers() + + assert np.allclose(P_base, P_new, rtol=1e-4) From e4ce7fb3e5846c933431d66feb9b1bb1419b0ccf Mon Sep 17 00:00:00 2001 From: misi9170 Date: Tue, 6 Aug 2024 10:36:18 -0600 Subject: [PATCH 31/37] Enable for all wake models/solvers. --- floris/core/solver.py | 70 ++++++++++++++++++++++ floris/utilities.py | 15 +++-- tests/data/input_full.yaml | 21 +++++++ tests/wind_direction_heterogeneity_test.py | 39 ++++++++++++ 4 files changed, 139 insertions(+), 6 deletions(-) diff --git a/floris/core/solver.py b/floris/core/solver.py index 8668e5c4f..2feded092 100644 --- a/floris/core/solver.py +++ b/floris/core/solver.py @@ -523,6 +523,20 @@ def cc_solver( z_i = np.mean(grid.z_sorted[:, i:i+1], axis=(2, 3)) z_i = z_i[:, :, None, None] + if grid.use_turbine_specific_layouts: + deficit_model_args = update_model_args( + deficit_model_args, + grid.x_sorted_per_turbine[:, :, :, :, i], + grid.y_sorted_per_turbine[:, :, :, :, i], + grid.z_sorted_per_turbine[:, :, :, :, i], + ) + deflection_model_args = update_model_args( + deflection_model_args, + grid.x_sorted_per_turbine[:, :, :, :, i], + grid.y_sorted_per_turbine[:, :, :, :, i], + grid.z_sorted_per_turbine[:, :, :, :, i], + ) + rotor_diameter_i = farm.rotor_diameters_sorted[:, i:i+1, None, None] mask2 = ( @@ -796,6 +810,20 @@ def full_flow_cc_solver( u_i = turbine_grid_flow_field.u_sorted[:, i:i+1] v_i = turbine_grid_flow_field.v_sorted[:, i:i+1] + if flow_field_grid.use_turbine_specific_layouts: + deficit_model_args = update_model_args( + deficit_model_args, + flow_field_grid.x_sorted_per_turbine[:, :, :, :, i], + flow_field_grid.y_sorted_per_turbine[:, :, :, :, i], + flow_field_grid.z_sorted_per_turbine[:, :, :, :, i], + ) + deflection_model_args = update_model_args( + deflection_model_args, + flow_field_grid.x_sorted_per_turbine[:, :, :, :, i], + flow_field_grid.y_sorted_per_turbine[:, :, :, :, i], + flow_field_grid.z_sorted_per_turbine[:, :, :, :, i], + ) + turb_avg_vels = average_velocity(turbine_grid_flow_field.u_sorted) turb_Cts = thrust_coefficient( velocities=turb_avg_vels, @@ -960,6 +988,20 @@ def turbopark_solver( z_i = np.mean(grid.z_sorted[:, i:i+1], axis=(2, 3)) z_i = z_i[:, :, None, None] + if grid.use_turbine_specific_layouts: + deficit_model_args = update_model_args( + deficit_model_args, + grid.x_sorted_per_turbine[:, :, :, :, i], + grid.y_sorted_per_turbine[:, :, :, :, i], + grid.z_sorted_per_turbine[:, :, :, :, i], + ) + deflection_model_args = update_model_args( + deflection_model_args, + grid.x_sorted_per_turbine[:, :, :, :, i], + grid.y_sorted_per_turbine[:, :, :, :, i], + grid.z_sorted_per_turbine[:, :, :, :, i], + ) + Cts = thrust_coefficient( velocities=flow_field.u_sorted, turbulence_intensities=flow_field.turbulence_intensity_field_sorted, @@ -1228,6 +1270,20 @@ def empirical_gauss_solver( z_i = np.mean(grid.z_sorted[:, i:i+1], axis=(2, 3)) z_i = z_i[:, :, None, None] + if grid.use_turbine_specific_layouts: + deficit_model_args = update_model_args( + deficit_model_args, + grid.x_sorted_per_turbine[:, :, :, :, i], + grid.y_sorted_per_turbine[:, :, :, :, i], + grid.z_sorted_per_turbine[:, :, :, :, i], + ) + deflection_model_args = update_model_args( + deflection_model_args, + grid.x_sorted_per_turbine[:, :, :, :, i], + grid.y_sorted_per_turbine[:, :, :, :, i], + grid.z_sorted_per_turbine[:, :, :, :, i], + ) + ct_i = thrust_coefficient( velocities=flow_field.u_sorted, turbulence_intensities=flow_field.turbulence_intensity_field_sorted, @@ -1452,6 +1508,20 @@ def full_flow_empirical_gauss_solver( z_i = np.mean(turbine_grid.z_sorted[:, i:i+1], axis=(2,3)) z_i = z_i[:, :, None, None] + if flow_field_grid.use_turbine_specific_layouts: + deficit_model_args = update_model_args( + deficit_model_args, + flow_field_grid.x_sorted_per_turbine[:, :, :, :, i], + flow_field_grid.y_sorted_per_turbine[:, :, :, :, i], + flow_field_grid.z_sorted_per_turbine[:, :, :, :, i], + ) + deflection_model_args = update_model_args( + deflection_model_args, + flow_field_grid.x_sorted_per_turbine[:, :, :, :, i], + flow_field_grid.y_sorted_per_turbine[:, :, :, :, i], + flow_field_grid.z_sorted_per_turbine[:, :, :, :, i], + ) + ct_i = thrust_coefficient( velocities=turbine_grid_flow_field.u_sorted, turbulence_intensities=turbine_grid_flow_field.turbulence_intensity_field_sorted, diff --git a/floris/utilities.py b/floris/utilities.py index bfc5b43a6..aafc086bf 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -619,14 +619,17 @@ def update_model_args( Update the x_sorted, y_sorted, and z_sorted entries of the model_args dict. """ # Check that new x, y, z are the same length as the old - if (x.shape != model_args["x"].shape or - y.shape != model_args["y"].shape or - z.shape != model_args["z"].shape): + if (("x" in model_args.keys() and x.shape != model_args["x"].shape) or \ + ("y" in model_args.keys() and y.shape != model_args["y"].shape) or \ + ("z" in model_args.keys() and z.shape != model_args["z"].shape)): raise ValueError("Size of new x,y,z arrays must match existing ones.") # Update the x, y, z values - model_args["x"] = x - model_args["y"] = y - model_args["z"] = z + if "x" in model_args.keys(): + model_args["x"] = x + if "y" in model_args.keys(): + model_args["y"] = y + if "z" in model_args.keys(): + model_args["z"] = z return model_args diff --git a/tests/data/input_full.yaml b/tests/data/input_full.yaml index 49c41273d..2ed76a01a 100644 --- a/tests/data/input_full.yaml +++ b/tests/data/input_full.yaml @@ -60,6 +60,12 @@ wake: ad: 0.0 bd: 0.0 kd: 0.05 + empirical_gauss: + horizontal_deflection_gain_D: 3.0 + vertical_deflection_gain_D: -1 + deflection_rate: 22 + mixing_gain_deflection: 0.0 + yaw_added_mixing_gain: 0.0 wake_velocity_parameters: cc: @@ -81,6 +87,19 @@ wake: turboparkgauss: A: 0.04 include_mirror_wake: True + turbopark: + A: 0.04 + sigma_max_rel: 4.0 + empirical_gauss: + wake_expansion_rates: + - 0.023 + - 0.008 + breakpoints_D: + - 10 + sigma_0_D: 0.28 + smoothing_length_D: 2.0 + mixing_gain_velocity: 2.0 + wake_turbulence_parameters: crespo_hernandez: @@ -88,3 +107,5 @@ wake: constant: 0.9 ai: 0.83 downstream: -0.25 + wake_induced_mixing: + atmospheric_ti_gain: 0.0 diff --git a/tests/wind_direction_heterogeneity_test.py b/tests/wind_direction_heterogeneity_test.py index 366493893..f23610138 100644 --- a/tests/wind_direction_heterogeneity_test.py +++ b/tests/wind_direction_heterogeneity_test.py @@ -205,3 +205,42 @@ def test_rotational_symmetry(): P_new = fmodel.get_turbine_powers() assert np.allclose(P_base, P_new, rtol=1e-4) + +def test_wake_models(): + fmodel = FlorisModel(configuration=YAML_INPUT) + fmodel_dict = fmodel.core.as_dict() + velocity_models = fmodel_dict["wake"]["wake_velocity_parameters"].keys() + + # Switch off secondary effects + fmodel_dict["wake"]["enable_secondary_steering"] = False + fmodel_dict["wake"]["enable_yaw_added_recovery"] = False + fmodel_dict["wake"]["enable_transverse_velocities"] = False + + # Loop through velocity models + for m in velocity_models: + fmodel_dict["wake"]["model_strings"]["velocity_model"] = m + if m == "empirical_gauss": + fmodel_dict["wake"]["model_strings"]["turbulence_model"] = "wake_induced_mixing" + fmodel_dict["wake"]["model_strings"]["deflection_model"] = "empirical_gauss" + elif m == "jensen": + fmodel_dict["wake"]["model_strings"]["turbulence_model"] = "crespo_hernandez" + fmodel_dict["wake"]["model_strings"]["deflection_model"] = "jimenez" + else: + fmodel_dict["wake"]["model_strings"]["turbulence_model"] = "crespo_hernandez" + fmodel_dict["wake"]["model_strings"]["deflection_model"] = "gauss" + + fmodel = FlorisModel(fmodel_dict) + + # First run without heterogeneity + fmodel.set(layout_x=layout_x, layout_y=layout_y) + fmodel.run() + P_without_het = fmodel.get_turbine_powers()[0,:] + + # Add wind direction heterogeneity and rerun + fmodel.set(heterogeneous_inflow_config=heterogeneous_wd_inflow_config_2D) + fmodel.run() + P_with_het = fmodel.get_turbine_powers()[0,:] + + # Confirm upstream power the same; downstream power increased + assert np.allclose(P_without_het[0], P_with_het[0]) + assert P_with_het[1] > P_without_het[1] From 86884bc789430c4a68b449d464f3dc2fdbeabf0d Mon Sep 17 00:00:00 2001 From: misi9170 Date: Tue, 6 Aug 2024 10:36:47 -0600 Subject: [PATCH 32/37] Rename to unit test. --- ...ogeneity_test.py => wind_direction_heterogeneity_unit_test.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename tests/{wind_direction_heterogeneity_test.py => wind_direction_heterogeneity_unit_test.py} (100%) diff --git a/tests/wind_direction_heterogeneity_test.py b/tests/wind_direction_heterogeneity_unit_test.py similarity index 100% rename from tests/wind_direction_heterogeneity_test.py rename to tests/wind_direction_heterogeneity_unit_test.py From 7f9daa8a81c4b9aa17aaf97b2660632d99ff3282 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Tue, 6 Aug 2024 14:34:34 -0600 Subject: [PATCH 33/37] Enable specifying wind speed heterogeneity using u and v --- floris/core/flow_field.py | 25 ++++++--- .../wind_direction_heterogeneity_unit_test.py | 52 ++++++++++++++++++- 2 files changed, 68 insertions(+), 9 deletions(-) diff --git a/floris/core/flow_field.py b/floris/core/flow_field.py index ce2d6529b..fc9fb4cc1 100644 --- a/floris/core/flow_field.py +++ b/floris/core/flow_field.py @@ -105,12 +105,18 @@ def heterogeneous_config_validator(self, instance: attrs.Attribute, value: dict return # Check that the correct keys are supplied for the heterogeneous_inflow_config dict - for k in ["speed_multipliers", "x", "y"]: + for k in ["x", "y"]: if k not in value.keys(): raise ValueError( - "heterogeneous_inflow_config must contain entries for 'speed_multipliers'," - f"'x', and 'y', with 'z' optional. Missing '{k}'." + "heterogeneous_inflow_config must contain entries for " + f"'x' and 'y', with 'z' optional. Missing '{k}'." ) + if ("speed_multipliers" not in value.keys() and + ("u" not in value.keys() or "v" not in value.keys())): + raise ValueError( + "heterogeneous_inflow_config must contain entries for either 'speed_multipliers'" + " or 'u' and 'v'." + ) if "z" not in value: # If only a 2D case, add "None" for the z locations value["z"] = None @@ -281,7 +287,13 @@ def generate_heterogeneous_wind_map(self): - **y**: A list of y locations at which the speed up factors are defined. - **z** (optional): A list of z locations at which the speed up factors are defined. """ - speed_multipliers = self.heterogeneous_inflow_config['speed_multipliers'] + if "speed_multipliers" in self.heterogeneous_inflow_config: + speed_multipliers = self.heterogeneous_inflow_config['speed_multipliers'] + else: + speed_multipliers = np.sqrt( + np.array(self.heterogeneous_inflow_config['u'])**2 + + np.array(self.heterogeneous_inflow_config['v'])**2 + ) / self.wind_speeds x = self.heterogeneous_inflow_config['x'] y = self.heterogeneous_inflow_config['y'] z = self.heterogeneous_inflow_config['z'] @@ -307,10 +319,7 @@ def generate_heterogeneous_wind_map(self): @property def _speed_heterogeneity(self): - if self.heterogeneous_inflow_config is not None: - return "speed_multipliers" in self.heterogeneous_inflow_config - else: - return False + return self.heterogeneous_inflow_config is not None @staticmethod def interpolate_multiplier_xy(x: NDArrayFloat, diff --git a/tests/wind_direction_heterogeneity_unit_test.py b/tests/wind_direction_heterogeneity_unit_test.py index f23610138..533e4e384 100644 --- a/tests/wind_direction_heterogeneity_unit_test.py +++ b/tests/wind_direction_heterogeneity_unit_test.py @@ -18,7 +18,7 @@ "y": [-1000, 1000, -1000, 1000], "u": [[8.0, 8.0, 8.0, 8.0]], "v": [[0.0, 0.0, 4.0, 4.0]], - "speed_multipliers": [[1.0, 1.0, 1.0, 1.0]] # Currently a necessary input + "speed_multipliers": [[1.0, 1.0, 1.0, 1.0]] # Removes wind speed variations } def test_power(): @@ -244,3 +244,53 @@ def test_wake_models(): # Confirm upstream power the same; downstream power increased assert np.allclose(P_without_het[0], P_with_het[0]) assert P_with_het[1] > P_without_het[1] + +def test_speed_multipliers_option(): + # Check that setting speed heterogeneity by u, v, works as expected + fmodel = FlorisModel(configuration=YAML_INPUT) + fmodel.set( + heterogeneous_inflow_config={ + "x": [0, 0, 1, 1], + "y": [0, 1, 0, 1], + "u": np.ones((1, 4)), + "v": np.ones((1, 4)), + }, + wind_shear=0.0 + ) + assert (fmodel.core.flow_field.u_initial_sorted == np.sqrt(2)).all() + + # Look at more realistic case + heterogeneous_wd_inflow_config_test = { + "x": [-1000, -1000, 1000, 1000], + "y": [-1000, 1000, -1000, 1000], + "u": [[8.0, 8.0, 8.0, 8.0]], + "v": [[0.0, 0.0, 4.0, 4.0]], + "speed_multipliers": [[1.0, 1.0, 1.0, 1.0]] + } + + # First run with speed_multipliers specified + fmodel.set( + layout_x=layout_x, + layout_y=layout_y, + heterogeneous_inflow_config=heterogeneous_wd_inflow_config_test + ) + fmodel.run() + P_with_sm = fmodel.get_turbine_powers()[0,:] + + # Same speed at second turbine + assert (fmodel.core.flow_field.u_initial_sorted[0,1,:,:] + == fmodel.core.flow_field.u_initial_sorted[0,0,:,:]).all() + + # Set without specified speed multipliers and rerun + # u, v will result in higher speeds at the back turbine + del heterogeneous_wd_inflow_config_test["speed_multipliers"] + fmodel.set(heterogeneous_inflow_config=heterogeneous_wd_inflow_config_test) + fmodel.run() + P_without_sm = fmodel.get_turbine_powers()[0,:] + + # Higher speeds at second turbine + assert (fmodel.core.flow_field.u_initial_sorted[0,1,:,:] + > fmodel.core.flow_field.u_initial_sorted[0,0,:,:]).all() + + # Also higher power than case without + assert P_without_sm[1] > P_with_sm[1] From e391664e78258671c27abc31b7fac0f8a58d9c2a Mon Sep 17 00:00:00 2001 From: misha Date: Thu, 8 Aug 2024 13:15:50 -0600 Subject: [PATCH 34/37] improved integration, handling for moving grid poitns by streamline. --- ...d_direction_heterogeneity_visualization.py | 2 +- floris/floris_model.py | 4 +- floris/utilities.py | 129 ++++++++++++++---- 3 files changed, 105 insertions(+), 30 deletions(-) diff --git a/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py index fb0d61b6d..54663169b 100644 --- a/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py +++ b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py @@ -13,7 +13,7 @@ fmodel = FlorisModel("../inputs/gch.yaml") -viz_by_sweep = True +viz_by_sweep = False # # Some wake models may not yet have a visualization method included, for these cases can use # # a slower version which scans a turbine model to produce the horizontal flow diff --git a/floris/floris_model.py b/floris/floris_model.py index 233ee0ba6..27b60649c 100644 --- a/floris/floris_model.py +++ b/floris/floris_model.py @@ -1041,15 +1041,13 @@ def set_for_viz(self, findex: int, solver_settings: dict) -> None: heterogeneous_inflow_config = { 'x': self.core.flow_field.heterogeneous_inflow_config['x'], 'y': self.core.flow_field.heterogeneous_inflow_config['y'], - 'speed_multipliers': - self.core.flow_field.heterogeneous_inflow_config['speed_multipliers'][findex:findex+1], } for k in ["z"]: # Optional, independent of findex if k in self.core.flow_field.heterogeneous_inflow_config: heterogeneous_inflow_config[k] = ( self.core.flow_field.heterogeneous_inflow_config[k] ) - for k in ["u", "v", "w"]: # Optional, dependent on findex + for k in ["u", "v", "speed_multipliers"]: # Optional, dependent on findex if k in self.core.flow_field.heterogeneous_inflow_config: if self.core.flow_field.heterogeneous_inflow_config[k] is not None: heterogeneous_inflow_config[k] = ( diff --git a/floris/utilities.py b/floris/utilities.py index aafc086bf..2f1a6e7a9 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -11,10 +11,11 @@ Tuple, ) +import numexpr as ne import numpy as np import yaml from attrs import define, field -from scipy.integrate import odeint +from scipy.integrate import solve_ivp from scipy.interpolate import LinearNDInterpolator from floris.type_dec import floris_array_converter, NDArrayFloat @@ -376,6 +377,7 @@ def apply_wind_direction_heterogeneity_warping( # TODO: Can I avoid any of these looping operations? Do I need to? for findex in range(n_findex): for tindex in range(n_turbines): + print("computing streamline for tindex", tindex) # 1. Compute the streamline streamline_start = np.array([ x_turbines[findex,tindex], @@ -390,8 +392,11 @@ def apply_wind_direction_heterogeneity_warping( ) # 2. Compute the point movements - x_shifted, y_shifted, z_shifted = shift_points_by_streamline( - streamline, x_points[findex,:,:,:], y_points[findex,:,:,:], z_points[findex,:,:,:] + x_shifted, y_shifted, z_shifted = shift_points_by_streamline2( + streamline, + x_points[findex,:,:,:], + y_points[findex,:,:,:], + z_points[findex,:,:,:], ) # 3. Apply to per_turbine values @@ -418,43 +423,55 @@ def compute_2D_streamline(x, y, u, v, xy_0, n_points=1000): uv = np.array([u,v]).T # Normalize - xy = (xy - offset) / scale - xy_0 = (xy_0 - offset) / scale + xy_normalized = (xy - offset) / scale + xy_0_normalized = (xy_0 - offset) / scale # TODO: Loop over findices! Maybe not necessary---can I compute all streamlines at once?? # Would it help? Currently assuming a single findex. - interp = LinearNDInterpolator(xy, uv) + interp = LinearNDInterpolator(xy_normalized, uv, fill_value=0) - def velocity_field_interpolate_reg(xy_, s): - if (xy_ < -0.5).any() or (xy_ > 0.5).any(): - return np.array([0, 0]) - else: - return interp(xy_)[0] + def velocity_field_interpolate_reg(s, xy_): + # if (xy_ < -0.5).any() or (xy_ > 0.5).any(): + # return np.array([0.1, 0.1])*np.sign(xy_) + # else: + return interp(xy_)[0] - streamline = odeint(velocity_field_interpolate_reg, xy_0, s) # n_points x 2 + ode_solution = solve_ivp( + velocity_field_interpolate_reg, + (0, 1), + xy_0_normalized, + t_eval=np.linspace(0,1,n_points) + ) + # Handle the ODE solution + streamline = ode_solution.y.T # Determine point of exit from specified domain if (streamline < -0.5).any(): low_exit = np.argwhere(streamline < -0.5)[0,0] else: - low_exit = 1000 + low_exit = n_points if (streamline > 0.5).any(): high_exit = np.argwhere(streamline > 0.5)[0,0] else: - high_exit = 1000 + high_exit = n_points exit_index = np.minimum(low_exit, high_exit) if exit_index == 0: raise ValueError("Streamline could not be propagated inside the domain.") - if exit_index == 1000: - raise ValueError("Streamline did not reach the edge of the domain.") + if exit_index == n_points: + # TODO: recommend ensuring domain covers the buffer region + if (abs(streamline) > 0.8*0.5).any(): # 10% buffer given, + pass + else: + import ipdb; ipdb.set_trace() + raise ValueError("Streamline did not reach the edge of the domain.") # Copy points after exit - streamline[exit_index:] = streamline[exit_index-1] + #streamline[exit_index:] = streamline[exit_index-1] # Denormalize - streamline = streamline * scale + offset + streamline_denormalized = streamline * scale + offset - return streamline + return streamline_denormalized def shift_points_by_streamline(streamline, x, y, z): ## I think should be dist by turbine only (not each rotor point); and then @@ -475,26 +492,86 @@ def shift_points_by_streamline(streamline, x, y, z): # TODO: for viz/solve for points, there are many more points than turbines. # It would be nice not to have a loop here. - for tindex in range(len(min_idx)): - streamline_diffs = np.diff(streamline[:min_idx[tindex], :], axis=0) + for i in range(len(min_idx)): + streamline_diffs = np.diff(streamline[:min_idx[i], :], axis=0) along_stream_dist = np.sum(np.linalg.norm(streamline_diffs, axis=1)) - off_stream_dist = dists_to_streamline[min_idx[tindex], tindex] - off_stream_disp = disps_to_streamline[min_idx[tindex], :, tindex] + off_stream_dist = dists_to_streamline[min_idx[i], i] + off_stream_disp = disps_to_streamline[min_idx[i], :, i] disp_x = off_stream_disp[0] disp_y = off_stream_disp[1] theta = np.arctan2(disp_y, disp_x) - x_new[tindex,:,:] = streamline[0,0] + along_stream_dist + x_new[i,:,:] = streamline[0,0] + along_stream_dist # y_new[tindex,:,:] = streamline[0,1] + disp_x*np.sin(theta) + disp_y*np.cos(theta)\ # + rotor_y[tindex,:,:] - y_new[tindex,:,:] = streamline[0,1] + np.sign(theta)*off_stream_dist + rotor_y[tindex,:,:] + y_new[i,:,:] = streamline[0,1] + np.sign(theta)*off_stream_dist + rotor_y[i,:,:] #y_new[tindex,:,:] = streamline[0,1] + off_stream_dist + rotor_y[tindex,:,:] #if tindex == 1: - # import ipdb; ipdb.set_trace() + return x_new, y_new, z + + #return x_new.reshape(x.shape), y_new.reshape(y.shape), z + +def shift_points_by_streamline2(streamline, x, y, z, rotor_points_only=True): + ## I think should be dist by turbine only (not each rotor point); and then + # shift all rotor points by the same amount. + # TODO: How will this work over findices? Need to work that out. + + #xy = np.concatenate((x.mean(axis=(1,2))[None,:],y.mean(axis=(1,2))[None,:]), axis=0) + xy = np.concatenate((x.flatten()[None,:], y.flatten()[None,:]), axis=0) + + #rotor_y = y - y.mean(axis=(1,2), keepdims=True) # TODO: still ok for viz solve? + + # Storage for new points + # x_new = x.copy() + # y_new = y.copy() + + disps_to_streamline = xy[None,:,:] - streamline[:, :,None] + # TODO: This line slow. What if norm is delayed? + dists_to_streamline = np.linalg.norm(disps_to_streamline, axis=1) + disps_x = disps_to_streamline[:, 0, :] + disps_y = disps_to_streamline[:, 1, :] + dists_to_streamline = ne.evaluate("sqrt(disps_x ** 2 + disps_y ** 2)") + min_idx = np.argmin(dists_to_streamline, axis=0) + #points = streamline[min_idx, :] + + # TODO: for viz/solve for points, there are many more points than turbines. + # It would be nice not to have a loop here. + # x_new = np.zeros(len(min_idx)) + # y_new = np.zeros(len(min_idx)) + + s_diffs = np.diff(streamline, axis=0) + as_dist_all = np.insert(np.cumsum(np.linalg.norm(s_diffs, axis=1)), 0, 0) + as_dists = as_dist_all[min_idx] + + disps = disps_to_streamline[min_idx, :, np.arange(0, xy.shape[1], 1)] + theta = np.arctan2(disps[:,1], disps[:,0]) + os_dists = dists_to_streamline[min_idx, np.arange(0, xy.shape[1], 1)] + # for i in range(len(min_idx)): + # streamline_diffs = np.diff(streamline[:min_idx[i], :], axis=0) + # along_stream_dist = np.sum(np.linalg.norm(streamline_diffs, axis=1)) + # off_stream_dist = dists_to_streamline[min_idx[i], i] + # off_stream_disp = disps_to_streamline[min_idx[i], :, i] + # disp_x = off_stream_disp[0] + # disp_y = off_stream_disp[1] + + # theta = np.arctan2(disp_y, disp_x) + + # x_new[i] = streamline[0,0] + along_stream_dist + # # y_new[tindex,:,:] = streamline[0,1] + disp_x*np.sin(theta) + disp_y*np.cos(theta)\ + # # + rotor_y[tindex,:,:] + # y_new[i] = streamline[0,1] + np.sign(theta)*off_stream_dist + # #y_new[tindex,:,:] = streamline[0,1] + off_stream_dist + rotor_y[tindex,:,:] + # #if tindex == 1: + # # import ipdb; ipdb.set_trace() + + x_new = (streamline[0,0] + as_dists).reshape(x.shape) + y_new = (streamline[0,1] + np.sign(theta)*os_dists).reshape(y.shape) return x_new, y_new, z + #return x_new.reshape(x.shape), y_new.reshape(y.shape), z + class Loader(yaml.SafeLoader): def __init__(self, stream): From 625c934868653060b82193f2102eeb4950fb3502 Mon Sep 17 00:00:00 2001 From: misha Date: Mon, 26 Aug 2024 10:38:40 -0600 Subject: [PATCH 35/37] Delay sqrt (doesnt seem to help much) --- floris/utilities.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/floris/utilities.py b/floris/utilities.py index 2f1a6e7a9..e76a7c410 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -531,8 +531,10 @@ def shift_points_by_streamline2(streamline, x, y, z, rotor_points_only=True): dists_to_streamline = np.linalg.norm(disps_to_streamline, axis=1) disps_x = disps_to_streamline[:, 0, :] disps_y = disps_to_streamline[:, 1, :] - dists_to_streamline = ne.evaluate("sqrt(disps_x ** 2 + disps_y ** 2)") - min_idx = np.argmin(dists_to_streamline, axis=0) + #dists_to_streamline = ne.evaluate("sqrt(disps_x ** 2 + disps_y ** 2)") + #min_idx = np.argmin(dists_to_streamline, axis=0) + dists_to_streamline_sq = disps_x ** 2 + disps_y ** 2 + min_idx = np.argmin(dists_to_streamline_sq, axis=0) #points = streamline[min_idx, :] # TODO: for viz/solve for points, there are many more points than turbines. @@ -546,7 +548,8 @@ def shift_points_by_streamline2(streamline, x, y, z, rotor_points_only=True): disps = disps_to_streamline[min_idx, :, np.arange(0, xy.shape[1], 1)] theta = np.arctan2(disps[:,1], disps[:,0]) - os_dists = dists_to_streamline[min_idx, np.arange(0, xy.shape[1], 1)] + #os_dists = dists_to_streamline[min_idx, np.arange(0, xy.shape[1], 1)] + os_dists = np.sqrt(dists_to_streamline_sq[min_idx, np.arange(0, xy.shape[1], 1)]) # for i in range(len(min_idx)): # streamline_diffs = np.diff(streamline[:min_idx[i], :], axis=0) # along_stream_dist = np.sum(np.linalg.norm(streamline_diffs, axis=1)) From e42b4f4315a89ed6a7f8b431a3ae5bc95a71db00 Mon Sep 17 00:00:00 2001 From: misi9170 Date: Thu, 15 Aug 2024 09:24:53 -0600 Subject: [PATCH 36/37] Switch viz_by_sweep to False. --- .../00x_wind_direction_heterogeneity_visualization.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py index 54663169b..f9a521bb0 100644 --- a/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py +++ b/examples/examples_heterogeneous/00x_wind_direction_heterogeneity_visualization.py @@ -18,7 +18,7 @@ # # Some wake models may not yet have a visualization method included, for these cases can use # # a slower version which scans a turbine model to produce the horizontal flow x = np.array([-100, 0, 100, 200, 300, 400, 500, 600]) -v = np.array([[0, 1, 2, 3, 4, 5, 6, 7]]) +v = -np.array([[0, 1, 2, 3, 4, 5, 6, 7]]) het_inflow_config = { "x": np.repeat(x, 2), From 607f50a16c7a50210d461107f65bd7be4f4c8a9e Mon Sep 17 00:00:00 2001 From: misi9170 Date: Mon, 26 Aug 2024 10:52:35 -0600 Subject: [PATCH 37/37] Revert to shift_points_by_streamline while debugging. --- floris/utilities.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/floris/utilities.py b/floris/utilities.py index e76a7c410..ef4e966e3 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -392,7 +392,7 @@ def apply_wind_direction_heterogeneity_warping( ) # 2. Compute the point movements - x_shifted, y_shifted, z_shifted = shift_points_by_streamline2( + x_shifted, y_shifted, z_shifted = shift_points_by_streamline( streamline, x_points[findex,:,:,:], y_points[findex,:,:,:],