diff --git a/examples/examples_curled_wake_model/01_test_cwm_solve.py b/examples/examples_curled_wake_model/01_test_cwm_solve.py new file mode 100644 index 000000000..d73ae5e5c --- /dev/null +++ b/examples/examples_curled_wake_model/01_test_cwm_solve.py @@ -0,0 +1,39 @@ +import numpy as np + +from floris import FlorisModel + + +fmodel = FlorisModel("../inputs/cwm.yaml") + +# Changing the wind farm layout uses FLORIS' set method to a two-turbine layout +fmodel.set(layout_x=[100, 500.0, 1000.0], layout_y=[0.0, 0.0, 0.0]) + + +single_condition = False + +if single_condition: + fmodel.set( + wind_directions=np.array([270.0]), + wind_speeds=np.array([8.0]), + turbulence_intensities=np.array([0.06]) + ) +else: + fmodel.set( + wind_directions=np.array([270.0, 180.0]), + wind_speeds=np.array([9.0, 8.0]), + turbulence_intensities=np.array([0.06, 0.06]) + ) + +fmodel.run() + +# There are functions to get either the power of each turbine, or the farm power +turbine_powers = fmodel.get_turbine_powers() / 1000.0 +farm_power = fmodel.get_farm_power() / 1000.0 + +print("Turbine powers") +print(turbine_powers) +print("Shape: ", turbine_powers.shape) + +print("Farm power") +print(farm_power) +print("Shape: ", farm_power.shape) diff --git a/examples/examples_curled_wake_model/02_random_layout.py b/examples/examples_curled_wake_model/02_random_layout.py new file mode 100644 index 000000000..1b9c99114 --- /dev/null +++ b/examples/examples_curled_wake_model/02_random_layout.py @@ -0,0 +1,93 @@ +import numpy as np + +from floris import FlorisModel + + +def main(): + # The FlorisModel class is the entry point for most usage. + # Initialize using an input yaml file + fmodel = FlorisModel("../inputs/cwm.yaml") + np.random.seed(0) + + # Changing the wind farm layout uses FLORIS' set method to a two-turbine layout + #fmodel.set(layout_x=[100, 500.0], layout_y=[300.0, 100.0]) + +# layout_x = [100, 100, 500, 600, 800, 1200] +# layout_y = [0, 100, 500, 0, 300, 400] + + # Example + D = 126. # Turbine diameter [m] + + # Test a random layout generator + N= 20 # Number of turbines + xd = 5 * D # Minimum distance in x-direction + yd = 3 * D # Minimum distance in y-direction + xlim = (0, np.sqrt(N*2.5) * xd) # X limits + ylim = (0, np.sqrt(N*2.5) * yd) # Y limits + layout_x, layout_y = random_layout(N=N, xd=5*D, yd=3*D, xlim=xlim, ylim=ylim) + print(layout_x, layout_y) + + yaw_angles = np.random.default_rng(0).uniform(-20, 20, size=(1, N)) + #yaw_angles = -20 * np.ones((1, N)) + print("Yaw angles:", yaw_angles) + + fmodel.set(layout_x=layout_x, layout_y=layout_y, + yaw_angles=yaw_angles, + ) + + + # Single wind condition + fmodel.set( + wind_directions=np.array([270.0]), + wind_speeds=np.array([8.0]), + turbulence_intensities=np.array([0.06]) + ) + + # After the set method, the run method is called to perform the simulation + fmodel.run() + + # There are functions to get either the power of each turbine, or the farm power + turbine_powers = fmodel.get_turbine_powers() / 1000.0 + farm_power = fmodel.get_farm_power() / 1000.0 + + print("Turbine powers") + print(turbine_powers) + print("Shape: ", turbine_powers.shape) + + print("Farm power") + print(farm_power) + print("Shape: ", farm_power.shape) + + + +def random_layout(N, xd, yd, xlim=(0, 100), ylim=(0, 100), max_attempts=100000): + layout_x = [] + layout_y = [] + attempts = 0 + + while len(layout_x) < N and attempts < max_attempts: + x_new = np.random.uniform(*xlim) + y_new = np.random.uniform(*ylim) + + # Check distances with existing points + if layout_x: + dx = np.abs(np.array(layout_x) - x_new) + dy = np.abs(np.array(layout_y) - y_new) + + # Valid if no other point is too close in BOTH directions + if np.any((dx < xd) & (dy < yd)): + attempts += 1 + continue + + layout_x.append(x_new) + layout_y.append(y_new) + attempts = 0 # reset since we succeeded + + if len(layout_x) < N: + raise RuntimeError(f"Could not place {N} points after {max_attempts} attempts.") + + return np.array(layout_x), np.array(layout_y) + + +if __name__ == "__main__": + main() diff --git a/examples/inputs/cwm.yaml b/examples/inputs/cwm.yaml new file mode 100644 index 000000000..89f777886 --- /dev/null +++ b/examples/inputs/cwm.yaml @@ -0,0 +1,66 @@ + +name: CWM +description: Three turbines using Curled Wake Model +floris_version: v4 + +logging: + console: + enable: true + level: WARNING + file: + enable: false + level: WARNING + +# TODO: Likely will need to overwrite this. +solver: + type: turbine_grid + turbine_grid_points: 3 + +farm: + layout_x: + - 0.0 + - 630.0 + - 1260.0 + layout_y: + - 0.0 + - 0.0 + - 0.0 + turbine_type: + - nrel_5MW + +flow_field: + air_density: 1.225 + reference_wind_height: -1 # -1 is code for use the hub height + turbulence_intensities: + - 0.06 + wind_directions: + - 270.0 + wind_shear: 0.12 + wind_speeds: + - 8.0 + wind_veer: 10.0 + +wake: + model_strings: + combination_model: none + deflection_model: none + turbulence_model: none + velocity_model: cwm + + enable_secondary_steering: false + enable_yaw_added_recovery: false + enable_transverse_velocities: false + enable_active_wake_mixing: false + + wake_deflection_parameters: + None: {} + + wake_velocity_parameters: + cwm: + C: 4.0 + dx_D: 0.08 + dy_D: 0.1 + dz_D: 0.1 + + wake_turbulence_parameters: + None: {} diff --git a/floris/core/__init__.py b/floris/core/__init__.py index e37f9c113..a292df393 100644 --- a/floris/core/__init__.py +++ b/floris/core/__init__.py @@ -47,6 +47,7 @@ from .wake import WakeModelManager from .solver import ( cc_solver, + curled_wake_solver, empirical_gauss_solver, full_flow_cc_solver, full_flow_empirical_gauss_solver, diff --git a/floris/core/core.py b/floris/core/core.py index 857b90fa4..73665b4de 100644 --- a/floris/core/core.py +++ b/floris/core/core.py @@ -12,6 +12,7 @@ from floris.core import ( BaseClass, cc_solver, + curled_wake_solver, empirical_gauss_solver, Farm, FlowField, @@ -198,6 +199,13 @@ def steady_state_atmospheric_condition(self): self.grid, self.wake ) + elif vel_model=="cwm": + curled_wake_solver( + self.farm, + self.flow_field, + self.grid, + self.wake + ) else: sequential_solver( self.farm, diff --git a/floris/core/solver.py b/floris/core/solver.py index 4325dc2ac..3b81618a3 100644 --- a/floris/core/solver.py +++ b/floris/core/solver.py @@ -24,9 +24,15 @@ yaw_added_turbulence_mixing, ) from floris.core.wake_velocity.empirical_gauss import awc_added_wake_mixing -from floris.type_dec import NDArrayFloat +from floris.type_dec import ( + floris_float_type, + NDArrayFloat, +) from floris.utilities import cosd +# For the CWM +from .solver_cwm import curled_wake_solver + def calculate_area_overlap(wake_velocities, freestream_velocities, y_ngrid, z_ngrid): """ diff --git a/floris/core/solver_cwm.py b/floris/core/solver_cwm.py new file mode 100644 index 000000000..7fdbbba29 --- /dev/null +++ b/floris/core/solver_cwm.py @@ -0,0 +1,266 @@ +import numpy as np +from scipy.ndimage.filters import gaussian_filter + +from floris.core import ( + Farm, + FlowField, + thrust_coefficient, + TurbineGrid, +) +from floris.core.wake import WakeModelManager +from floris.type_dec import ( + floris_float_type, + NDArrayFloat, +) + + +def curled_wake_solver( + farm: Farm, + flow_field: FlowField, + grid: TurbineGrid, + model_manager: WakeModelManager +) -> NDArrayFloat: + # In the normal paradigm, flow_field contains quantities over the turbine + # rotors (perhaps not the best name, but that's how it is). + + # I propose that we start by keeping it that way; i.e., flow_field.u_sorted + # is the "final" output of the solve, with dimensions (n_findex, n_turbines, n_grid, n_grid), + # with n_grid being a number of grid points used for averaging + # (potentially urelated to the CWM grid). + + # I'm generating the data matrices for the CWM with an empty first dimension, + # as a placeholder for a final implementation with an n_findex dimension. + + dx = model_manager.velocity_model.dx_D * farm.rotor_diameters_sorted.min() + dy = model_manager.velocity_model.dy_D * farm.rotor_diameters_sorted.min() + dz = model_manager.velocity_model.dz_D * farm.rotor_diameters_sorted.min() + + #dx /= 2 # ensure smaller dx for numerical stability + x_min = 0 - farm.rotor_diameters_sorted.max() # 1D upstream of first turbine + x_max = grid.x_sorted.max() + farm.rotor_diameters_sorted.max() # 1D downstream of last turbine + y_min = grid.y_sorted.min() - 5*farm.rotor_diameters_sorted.max() + y_max = grid.y_sorted.max() + 5*farm.rotor_diameters_sorted.max() + x_1d = np.arange(x_min, x_max, dx, dtype=floris_float_type) # x coordinates + y_1d = np.arange(y_min, y_max, dy, dtype=floris_float_type) # y coordinates + z_1d = np.arange( + 0, farm.hub_heights_sorted.max()+ 2*farm.rotor_diameters_sorted.max(), + dz, + dtype=floris_float_type + ) # Go up to 0.5D above highest tip point + + n_x_planes = x_1d.shape[0] + + # Create large arrays (memory intensive!) + x, y, z = np.meshgrid(x_1d, y_1d, z_1d, indexing='ij') + x = np.repeat(x[None,:,:,:], flow_field.n_findex, axis=0) + y = np.repeat(y[None,:,:,:], flow_field.n_findex, axis=0) + z = np.repeat(z[None,:,:,:], flow_field.n_findex, axis=0) + + # Use the same number of x and y points for generality once we get to multiple + # findices at once (can consider revising later) + # Not yet handling heterogeneity; will work that in later. + # TODO: what is the 3 for? Does the solve not work for lower wind speeds? We will need to handle + # this better, if so. + u_freestream = ( + np.ones_like(x, dtype=floris_float_type) + * flow_field.wind_speeds[:, None, None, None] + ) + v_freestream = np.zeros_like(x, dtype=floris_float_type) # May not be needed. Always zeros? + w_freestream = np.zeros_like(x, dtype=floris_float_type) # May not be needed. Always zeros? + u_waked = np.zeros_like(x, dtype=floris_float_type) + v_waked = np.zeros_like(x, dtype=floris_float_type) + w_waked = np.zeros_like(x, dtype=floris_float_type) + + # + # Add flow features + # + + # Add boundary layer + print('Adding boundary layer:') + u_sheared = (z / flow_field.reference_wind_height) ** flow_field.wind_shear * u_freestream + u_freestream = np.maximum(u_sheared, u_freestream*.3) # Ensure freestream is above 3 m/s + + # Add veer + if np.abs(flow_field.wind_veer) > 1e-6: + print("Adding veer:") + v_sheared = model_manager.velocity_model.add_veer_theta( + z, + u_freestream, + theta_deg=flow_field.wind_veer, + th=farm.hub_heights_sorted.max(), + D=farm.rotor_diameters_sorted.max() + ) + print("V shear shape:", v_sheared.shape) + v_freestream += v_sheared + + # Compute the numerical viscosity needed for stability + Re = model_manager.velocity_model.Re + # Minimum viscosity for numerical stability + nu_min = u_freestream * flow_field.reference_wind_height / Re + # Initialize the turbulence viscosity based on the standard curled wake model formulation + nu_t = model_manager.velocity_model.initialize_turbulence_viscosity(u_freestream, z, dz, nu_min) + + # Code to generate one or multiple turbine indices that correspond to a certain x location + # Each element in the list turbines_in_plane should be a list of the turbines that appear at + # that plane's x location (i.e., often an empty list). + turbine_plane_map = np.zeros((flow_field.n_findex, farm.n_turbines), dtype=int) + + # Let's try to extract the turbines in each plane + for t in range(farm.n_turbines): + # Get the x coordinate of the turbine + turbine_x = np.mean(grid.x_sorted[:, t, :, :], axis=(1,2)) + # Find the index of the plane that contains this turbine + ip = np.argmin(np.abs(x_1d.reshape(-1,1) - turbine_x.reshape(1,-1)), axis=0) + turbine_plane_map[:, t] = ip + + count_inner_loops = 0 # Temporary + + for i in range(1, n_x_planes): + # Ensure freestream velocity is above 3 m/s (should really be 20% of U_inf) + u_freestream_plane = np.maximum(u_freestream[:, i, :, :], 3) # TODO: 3 m/s minimum + v_freestream_plane = v_freestream[:, i, :, :] + w_freestream_plane = w_freestream[:, i, :, :] + u_waked_plane = u_waked[:, i, :, :] + v_waked_plane = v_waked[:, i, :, :] # noqa: F841 TODO: remove if not used + w_waked_plane = w_waked[:, i, :, :] # noqa: F841 TODO: remove if not used + + # Extract the 2D slices for the current plane + u_prev = u_waked[:, i-1, :, :] # shape (n_findex, ny, nz) + v_prev = v_waked[:, i-1, :, :] # shape (n_findex, ny, nz) + w_prev = w_waked[:, i-1, :, :] # shape (n_findex, ny, nz) + u_fs = u_freestream_plane # shape (n_findex, ny, nz) + v_fs = v_freestream_plane # shape (n_findex, ny, nz) + w_fs = w_freestream_plane # shape (n_findex, ny, nz) + nu_2d = nu_t[:, i, :, :] # shape (n_findex, ny, nz) + + # Run Runge-Kutta step + u_new = model_manager.velocity_model.function( + u_prev, + dx, + u_fs, + v_fs + v_prev, + w_fs + w_prev, + dy, + dz, + nu_2d + ) + + # Assign it back + u_waked[:, i, :, :] = u_new + + # A simple evolution model to scale v the same way that U has scaled + # This saves all the work of resolving the transport equation (du/dx~dv/dx) + U = u_freestream_plane[0, :, :] + fact = (U + u_waked[0, i-1, :, :]) / (U + u_waked[0, i, :, :]) + # This ensures that V and W do not become larger (they should always decay) + fact = np.clip(fact, .1, 1.) + v_waked[:,i,:,:] = v_waked[:, i-1, :, :] * fact + w_waked[:,i,:,:] = w_waked[:, i-1, :, :] * fact + + # For this plane, find the findex-turbines that are present and modify the flow + findex_turbine_indices = np.argwhere(turbine_plane_map == i) + for fti in findex_turbine_indices: + count_inner_loops += 1 # Temporary + + f, t = fti # f is the findex, t is the turbine index + + # Let's get the rotor location for now (should access properly later) + t_y = np.mean(grid.y_sorted[f,t,:,:]) + t_z = np.mean(grid.z_sorted[f,t,:,:]) + rotor_diameter_t = farm.rotor_diameters_sorted[f, t] + + # Create a mask for points inside the rotor disk at this x-plane + # mask shape: (n_findex, n_y, n_z), True where (y, z) is inside rotor + rotor_mask = ( + ((y[f, i, :, :] - t_y) ** 2 + (z[f, i, :, :] - t_z) ** 2) + < (rotor_diameter_t / 2)**2 + ) + #print("Rotor mask shape:", rotor_mask.shape) + # The filtered mask for points inside the rotor disk (wider) + rotor_mask_filt = ( + ((y[f, i, :, :] - t_y) ** 2 + (z[f, i, :, :] - t_z) ** 2) + < (1.3 * rotor_diameter_t / 2)**2 + ) + # TODO: where does this 1.3 factor come from? Is it a user setting? + + u_rotor_values = u_freestream_plane[f, rotor_mask] + u_waked_plane[f, rotor_mask] + u_rotor_disk = np.mean(u_rotor_values) # TODO: Consider shear? + flow_field.u_sorted[f, t, :, :] = u_rotor_disk # or shape match if 2D needed + + ct = thrust_coefficient( + velocities=flow_field.u_sorted, + turbulence_intensities=flow_field.turbulence_intensity_field_sorted, + air_density=flow_field.air_density, + yaw_angles=farm.yaw_angles_sorted, + tilt_angles=farm.tilt_angles_sorted, + power_setpoints=farm.power_setpoints_sorted, + awc_modes=farm.awc_modes_sorted, + awc_amplitudes=farm.awc_amplitudes_sorted, + thrust_coefficient_functions=farm.turbine_thrust_coefficient_functions, + tilt_interps=farm.turbine_tilt_interps, + correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted, + turbine_type_map=farm.turbine_type_map_sorted, + turbine_power_thrust_tables=farm.turbine_power_thrust_tables, + ix_filter=[t], + average_method=grid.average_method, + cubature_weights=grid.cubature_weights, + multidim_condition=flow_field.multidim_conditions, + )[f] + + a = (1. - np.sqrt(1. - ct * np.cos(0)**2)) / 2 + + # TODO: is there a reason not to use the axial induction method? + # force scalar # Set a limit to guarantee numerical stability + a = float(np.minimum(a, 0.35)) # + # print("a:", a) + + # Apply induction to points inside the rotor disk + u_waked_plane[f, rotor_mask] = -2 * a * u_freestream_plane[f, rotor_mask] + #u_waked_plane = gaussian_filter(u_waked_plane, 3) + # (this filtered the entire plane) + # u_waked[0, i, :, :] = gaussian_filter(u_waked_plane, 3) + # Filter only around the rotor disk + u_waked[f,i,:,:][rotor_mask_filt] = gaussian_filter( + u_waked_plane[f, rotor_mask_filt], + 2 + ) + + # Add Curl + vcurl, wcurl = model_manager.velocity_model.add_curl( + y[:, i, :, :]-t_y, #[rotor_mask_filt], + z[:, i, :, :]-t_z, #[rotor_mask_filt], + D=rotor_diameter_t, + Ct=ct, + Uh=u_rotor_disk, + alpha=np.deg2rad(farm.yaw_angles_sorted[f, t]), + tilt=np.deg2rad(farm.tilt_angles_sorted[f, t]), + ground=True, + th=farm.hub_heights_sorted[f, t], + N=20, + ) + v_waked[:,i,:,:] += vcurl + w_waked[:,i,:,:] += wcurl + + # Ensure boundary conditions are satisfied + u_waked[:, i, :, [0, -1]] = 0 + u_waked[:, i, [0, -1], :] = 0 + + # This is the end of the loop that goes through the x planes. + # k = np.argmin(np.abs(z_1d - flow_field.reference_wind_height)) + #print("z_1d:", z_1d) + #print("dx:", dx) + # #print("k:", k) + # import matplotlib.pyplot as plt + + # #plt.pcolormesh(x_1d, y_1d, u_freestream[0, :, :, k].T, label='Waked', shading='gouraud') + # plt.pcolormesh(x_1d, y_1d, u_waked[0, :, :, k].T, label='Waked', shading='gouraud', + # #vmin=3, vmax=9, + # ) + # plt.colorbar() + # plt.gca().set_aspect('equal', adjustable='box') + # plt.show() + # print("Solve complete.") + + # Result: flow_field.u_sorted. + #print("Inner loop ran {0} times".format(count_inner_loops)) + return None diff --git a/floris/core/wake.py b/floris/core/wake.py index e58a85cb4..8de4ec387 100644 --- a/floris/core/wake.py +++ b/floris/core/wake.py @@ -6,6 +6,7 @@ from floris.core.wake_combination import ( FLS, MAX, + NoneWakeCombination, SOSFS, ) from floris.core.wake_deflection import ( @@ -21,6 +22,7 @@ ) from floris.core.wake_velocity import ( CumulativeGaussCurlVelocityDeficit, + CurledWakeVelocityDeficit, EmpiricalGaussVelocityDeficit, GaussVelocityDeficit, JensenVelocityDeficit, @@ -32,6 +34,7 @@ MODEL_MAP = { "combination_model": { + "none": NoneWakeCombination, "fls": FLS, "max": MAX, "sosfs": SOSFS @@ -50,6 +53,7 @@ "velocity_model": { "none": NoneVelocityDeficit, "cc": CumulativeGaussCurlVelocityDeficit, + "cwm": CurledWakeVelocityDeficit, "gauss": GaussVelocityDeficit, "jensen": JensenVelocityDeficit, "turbopark": TurbOParkVelocityDeficit, diff --git a/floris/core/wake_combination/__init__.py b/floris/core/wake_combination/__init__.py index 246aab65c..3deee0674 100644 --- a/floris/core/wake_combination/__init__.py +++ b/floris/core/wake_combination/__init__.py @@ -1,4 +1,5 @@ from floris.core.wake_combination.fls import FLS from floris.core.wake_combination.max import MAX +from floris.core.wake_combination.none import NoneWakeCombination from floris.core.wake_combination.sosfs import SOSFS diff --git a/floris/core/wake_combination/none.py b/floris/core/wake_combination/none.py new file mode 100644 index 000000000..a8abaeca8 --- /dev/null +++ b/floris/core/wake_combination/none.py @@ -0,0 +1,29 @@ + +from typing import Any, Dict + +import numpy as np +from attrs import define, field + +from floris.core import BaseModel + + +@define +class NoneWakeCombination(BaseModel): + """ + The None wake turbulence model is a placeholder code that simple ignores + any combination and returns None. + """ + + def prepare_function(self) -> dict: + pass + + def function( + self, + ambient_TI: float, + x: np.ndarray, + x_i: np.ndarray, + rotor_diameter: float, + axial_induction: np.ndarray, + ) -> None: + """Return None""" + return None diff --git a/floris/core/wake_velocity/__init__.py b/floris/core/wake_velocity/__init__.py index 07762a5cd..46c4b7364 100644 --- a/floris/core/wake_velocity/__init__.py +++ b/floris/core/wake_velocity/__init__.py @@ -1,5 +1,6 @@ from floris.core.wake_velocity.cumulative_gauss_curl import CumulativeGaussCurlVelocityDeficit +from floris.core.wake_velocity.curled_wake import CurledWakeVelocityDeficit from floris.core.wake_velocity.empirical_gauss import EmpiricalGaussVelocityDeficit from floris.core.wake_velocity.gauss import GaussVelocityDeficit from floris.core.wake_velocity.jensen import JensenVelocityDeficit diff --git a/floris/core/wake_velocity/curled_wake.py b/floris/core/wake_velocity/curled_wake.py new file mode 100644 index 000000000..2b2d748c3 --- /dev/null +++ b/floris/core/wake_velocity/curled_wake.py @@ -0,0 +1,309 @@ + +from typing import Any, Dict + +import numexpr as ne +import numpy as np +from attrs import ( + define, + field, + fields, +) + +from floris.core import ( + BaseModel, + Farm, + FlowField, + Grid, + Turbine, +) +from floris.type_dec import floris_float_type + + +NUM_EPS = fields(BaseModel).NUM_EPS.default + +@define +class CurledWakeVelocityDeficit(BaseModel): + """ + TO DO + """ + + C: float = field(default=4) # Turbulence constant for the CWM + dx_D: float = field(default=0.08) # grid size dx/D, recommended to be slightly smaller than dy/D + dy_D: float = field(default=0.1) # grid size dy/D, recommended to be around 0.1 + dz_D: float = field(default=0.1) # grid size dz/D, recommended to be the same as dy/D + Re: float = field(default=1e4) # Reynolds number for numerical stability + + def prepare_function( + self, + grid: Grid, + flow_field: FlowField, + ) -> Dict[str, Any]: + """ + This function prepares the inputs from the various FLORIS data structures + for use in the Jensen model. This should only be used to 'initialize' + the inputs. For any data that should be updated successively, + do not use this function and instead pass that data directly to + the model function. + """ + kwargs = { + "x": grid.x_sorted, + "y": grid.y_sorted, + "z": grid.z_sorted, + } + return kwargs + + # @profile + def function( + self, + u_prev: np.ndarray, + dx: float, + u_fs: np.ndarray, + v_fs: np.ndarray, + w_fs: np.ndarray, + dy: float, + dz: float, + nu_2d: np.ndarray, + ) -> None: + """ + Solve a plane of the curled wake model. + Return the velocity deficit at the next plane. + """ + + u_new = runge_kutta_step( + u_prev, + dx, + u_fs, + v_fs, + w_fs, + dy, + dz, + nu_2d + ) + + return u_new + + def initialize_turbulence_viscosity(self, u_freestream, Z, dz, nu_min, f=1.0, C=4.0): + """ + Initialize the turbulence viscosity based on the freestream velocity. + This is a simple initialization that can be modified as needed. + """ + dudz = np.gradient(u_freestream, dz, axis=-1) # Gradient in the z-direction + + # von-Karman constant + kappa = 0.41 + + # Mixing length based on: + # ALFRED K. BLACKADAR 1962 + # JIELUN SUN 2011 + # ~ lmda = f * 15. # The maximum length [m] + lmda = f * 27. # The maximum length [m] + lm = kappa * Z / (1 + kappa * Z / lmda) + + # Turbulent viscosity + nu_t = lm**2 * np.abs(dudz) + + # Pick the maximum between the turbulent model and the one required for + # numerical stability + nu = np.maximum(C * nu_t, nu_min) + + return nu + + def add_curl(self, + Y, Z, *, # <- no V, W here + D, Ct, Uh, alpha, tilt, + ground=False, th=0.0, + N=20, eps=None + ): + """ + Return ONLY the curl-induced velocities (Vcurl, Wcurl), without + modifying any inputs. Caller can add these to their own V, W. + """ + if eps is None: + eps = 0.2 * D + + # Lamb–Oseen-style kernel; drop-in consistent with your original. + def vortex(dy, dz, Gamma=1.0, eps=0.2): + r2 = dy*dy + dz*dz + # zero exactly at r=0; avoids masked assignment pitfalls + safe = r2 > 1e-12 + factor = np.zeros_like(r2, dtype=float) + np.divide(1.0 - np.exp(-r2/(eps*eps)), r2, out=factor, where=safe) + factor *= (float(Gamma) / (2.0*np.pi)) + uy = factor * dz # +dz + uz = -factor * dy # -dy + return uy, uz + + Vcurl = np.zeros_like(Y, dtype=float) + Wcurl = np.zeros_like(Z, dtype=float) + + R = 0.5 * D + + # Rotor normal + n = np.array([ + np.cos(alpha) * np.cos(tilt), # x + np.sin(alpha), # y + -np.cos(alpha) * np.sin(tilt) # z + ], dtype=float) + + inflow = np.array([1.0, 0.0, 0.0], dtype=float) + + # Total rotation angle + cosang = float(np.clip(np.dot(n, inflow), -1.0, 1.0)) + theta_total = np.arccos(cosang) + + # Rotation axis (use 1e-6 threshold like your original) + axis = np.cross(inflow, n) + axis_norm = np.linalg.norm(axis) + if axis_norm < 1e-6: + axis = np.array([0.0, 0.0, 1.0], dtype=float) + else: + axis /= axis_norm + axis_y, axis_z = axis[1], axis[2] + + # Circulation magnitude (same as original) + Gamma_total = -( + np.pi * D / 4.0 * 0.5 * Ct * Uh * + np.sin(theta_total) * np.cos(theta_total)**2 + ) + Gamma0 = 4.0 / np.pi * Gamma_total + + # Midpoint integration in θ + theta_edges = np.linspace(0.0, 0.5*np.pi, N + 1) + dtheta = theta_edges[1] - theta_edges[0] + theta = theta_edges[:-1] + 0.5*dtheta + r = R * np.sin(theta) + dr = R * np.cos(theta) * dtheta + + for s, ds in zip(r, dr): + # IMPORTANT: replicate original clipping order + denom = np.sqrt(1.0 - (2.0*s/D)**2) + denom = max(denom, 1e-12) + + Gamma = -4.0 * Gamma0 * s * ds / (D*D * denom) + + off_y, off_z = s * axis_y, s * axis_z + + # Primary pair + vt1, wt1 = vortex(Y - off_y, Z - off_z, Gamma, eps) + vt2, wt2 = vortex(Y + off_y, Z + off_z, -Gamma, eps) + Vcurl += vt1 + vt2 + Wcurl += wt1 + wt2 + + if ground: + z_offset = 2.0 * th + # Image vortices (signs match your original) + vt1g, wt1g = vortex(Y - off_y, Z + off_z + z_offset, -Gamma, eps) + vt2g, wt2g = vortex(Y + off_y, Z - off_z + z_offset, Gamma, eps) + Vcurl += vt1g + vt2g + Wcurl += wt1g + wt2g + + return Vcurl, Wcurl + + def add_veer_theta(self, Z, U, *, theta_deg, th, D): + """ + Compute the veer-induced spanwise velocity (theta-mode only). + + Parameters + ---------- + Z : array_like + Height array (same shape as U). + U : array_like or float + Streamwise velocity at each Z (broadcastable to Z). + (Use your hub or local profile; matches your original use of self.U.) + theta_deg : float + Total veer angle [degrees] between bottom and top of the rotor. + th : float + Hub height [m]. + D : float + Rotor diameter [m]. + + Returns + ------- + V_veer : ndarray + Spanwise velocity induced by veer (same shape as Z). + """ + # Linear angle profile centered at hub height: + # alpha(z) varies from +theta/2 at z = th - D/2 to -theta/2 at z = th + D/2, + # and extrapolates linearly outside that range. + # + # Closed-form equivalent of the interp1d in your code: + # angle_half = deg2rad(theta)/2 + # alpha(z) = - (deg2rad(theta) / D) * (z - th) + alpha = - np.deg2rad(theta_deg) * (Z - float(th)) / float(D) + + # Veer crossflow = tan(local angle) * local streamwise speed + V_veer = np.tan(alpha) * U + return V_veer + + + + + +def finite_diff(arr, dy, dz): + """ + Compute finite differences in a grid with 'ij' indexing: + - First dimension (axis 0): y-direction + - Second dimension (axis 1): z-direction + + Parameters: + arr: 2D numpy array of shape (ny, nz) + dy: Grid spacing in the y-direction (axis 0) + dz: Grid spacing in the z-direction (axis 1) + + Returns: + duwdy: Approximation of partial derivative w.r.t. y (axis 0) + duwdz: Approximation of partial derivative w.r.t. z (axis 1) + """ + duwdy = np.zeros_like(arr, dtype=floris_float_type) # Partial derivative w.r.t. y + duwdz = np.zeros_like(arr, dtype=floris_float_type) # Partial derivative w.r.t. z + + # Central difference (interior points) + duwdy[:,1:-1,:] = (arr[:,2:,:] - arr[:,:-2,:]) / (2 * dy) + duwdz[:,:,1:-1] = (arr[:,:,2:] - arr[:,:,:-2]) / (2 * dz) + + return duwdy, duwdz + + + +def laplacian(u, dy, dz): + lap = np.zeros_like(u, dtype=np.float64) + + # Compute second derivatives in y-direction (axis=1) + lap[:,2:-2,2:-2] = ( + (u[:,:-4,2:-2] - 2 * u[:,2:-2,2:-2] + u[:,4:,2:-2]) / (4 * dy * dy) # y-direction + + (u[:,2:-2,:-4] - 2 * u[:,2:-2,2:-2] + u[:,2:-2,4:]) / (4 * dz * dz) # z-direction + ) + + lap[:,0,:] = lap[:,0,:] + (u[:,2,:]/2 - u[:,0,:]/2 - u[:,1,:] + u[:,0,:]) / (dy * dy) + lap[:,1,:] = lap[:,1,:] + (u[:,3,:]/2 - u[:,1,:]/2 - u[:,1,:] + u[:,0,:]) / (2 * dy * dy) + lap[:,-2,:] = lap[:,-2,:] + (u[:,-1,:] - u[:,-2,:] - u[:,-2,:]/2 + u[:,-4,:]/2) / (2 * dy * dy) + lap[:,-1,:] = lap[:,-1,:] + (u[:,-1,:] - u[:,-2,:] - u[:,-1,:]/2 + u[:,-3,:]/2) / (dy * dy) + + lap[:,:,0] = lap[:,:,0] + (u[:,:,2]/2 - u[:,:,0]/2 - u[:,:,1] + u[:,:,0]) / (dz * dz) + lap[:,:,1] = lap[:,:,1] + (u[:,:,3]/2 - u[:,:,1]/2 - u[:,:,1] + u[:,:,0]) / (2 * dz * dz) + lap[:,:,-2] = lap[:,:,-2] + (u[:,:,-1] - u[:,:,-2] - u[:,:,-2]/2 + u[:,:,-4]/2) / (2 * dz * dz) + lap[:,:,-1] = lap[:,:,-1] + (u[:,:,-1] - u[:,:,-2] - u[:,:,-1]/2 + u[:,:,-3]/2) / (dz * dz) + + return lap + + +def compute_rhs_steady(u_current, U, V, W, dy, dz, nu): + duwdy, duwdz = finite_diff(u_current, dy, dz) + inv_U = 1.0 / (U + u_current) # Precompute inverse + rhs = inv_U * (-V * duwdy - W * duwdz + nu * laplacian(u_current, dy, dz)) + return rhs + + +def runge_kutta_step(u_current, dx, U, V, W, dy, dz, nu): + k1 = compute_rhs_steady(u_current, U, V, W, dy, dz, nu) + + tmp = u_current + 0.5 * dx * k1 # In-place variable reuse + k2 = compute_rhs_steady(tmp, U, V, W, dy, dz, nu) + + tmp = u_current + 0.5 * dx * k2 + k3 = compute_rhs_steady(tmp, U, V, W, dy, dz, nu) + + tmp = u_current + dx * k3 + k4 = compute_rhs_steady(tmp, U, V, W, dy, dz, nu) + + return u_current + (dx / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4) diff --git a/pyproject.toml b/pyproject.toml index 3c62ccc2b..359d052a2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -40,6 +40,7 @@ dependencies = [ "shapely~=2.0", "coloredlogs~=15.0", "pathos~=0.3", + "numba", ] [project.optional-dependencies]