Merge pull request #78 from byuflowlab/typestability

Typestability

Author: BTV25

src/fatigue_model.jl

@@ -645,9 +645,9 @@ function get_single_damage_surr(turbine_x,turbine_y,turbine_z,rotor_diameter,hub
 
         temp_resource = DiscretizedWindResource([3*pi/2], [ws], [1.0], measurementheight, air_density, ambient_tis, wind_shear_model)
 
-        turbine_velocities, turbine_ct, turbine_local_ti = turbine_velocities_one_direction(turbine_x, turbine_y, turbine_z, rotor_diameter, hub_height, turbine_yaw,
-                            turbine_ai, sorted_turbine_index, ct_model, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
-                            model_set; wind_farm_state_id=state_ID)
+        turbine_velocities, turbine_ct, _, turbine_local_ti = turbine_velocities_one_direction_full(turbine_x, turbine_y, turbine_z, rotor_diameter, hub_height, turbine_yaw,
+                            sorted_turbine_index, ct_model, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
+                            model_set; wind_farm_state_id=state_ID, turbine_ai=turbine_ai)
 
 
         for i = 1:nCycles*naz

src/general_models.jl

@@ -10,8 +10,7 @@ author: Benjamin Varela
 build_wind_farm_struct(x,turbine_x,turbine_y,turbine_z,hub_height,turbine_yaw,rotor_diameter,
             ct_models,generator_efficiency,cut_in_speed,cut_out_speed,rated_speed,rated_power,wind_resource,
             power_models,model_set,update_function;rotor_sample_points_y=[0.0],rotor_sample_points_z=[0.0],
-            AEP_scale=0.0,input_type=nothing,opt_x=false,opt_y=false,opt_hub=false,opt_yaw=false,opt_diam=false,
-            force_single_thread=false)
+            AEP_scale=0.0,input_type=nothing,opt_x=false,opt_y=false,opt_hub=false,opt_yaw=false,opt_diam=false)
 
 function to build a wind_farm_struct
 
@@ -42,17 +41,13 @@ function to build a wind_farm_struct
 - `opt_hub`: Boolean to optimize hub heights of turbines
 - `opt_yaw`: Boolean to optimize yaw angles of turbines
 - `opt_diam`: Boolean to optimize rotor diameters of turbines
-- `force_single_thread`: Boolean to force single thread calculation
 """
 function build_wind_farm_struct(x,turbine_x,turbine_y,turbine_z,hub_height,turbine_yaw,rotor_diameter,
             ct_models,generator_efficiency,cut_in_speed,cut_out_speed,rated_speed,rated_power,wind_resource,
             power_models,model_set,update_function;rotor_sample_points_y=[0.0],rotor_sample_points_z=[0.0],
-            AEP_scale=0.0,input_type=nothing,opt_x=false,opt_y=false,opt_hub=false,opt_yaw=false,opt_diam=false,
-            force_single_thread=false)
+            AEP_scale=0.0,input_type=nothing,opt_x=false,opt_y=false,opt_hub=false,opt_yaw=false,opt_diam=false)
 
     n_turbines = length(turbine_x)
-    n_threads = Threads.nthreads()
-    force_single_thread && (n_threads = 1)
     results = DiffResults.GradientResult(x)
     AEP_gradient = zeros(eltype(x),length(x))
     AEP = Array{eltype(x),0}(undef)
@@ -65,6 +60,7 @@ function build_wind_farm_struct(x,turbine_x,turbine_y,turbine_z,hub_height,turbi
                 cut_out_speed, rated_speed, rated_power, wind_resource, power_models, model_set;
                 rotor_sample_points_y=rotor_sample_points_y, rotor_sample_points_z=rotor_sample_points_z)
 
+    ideal_AEP == 0.0 && (ideal_AEP = 1.0)
     (AEP_scale == 0.0) && (AEP_scale = 1.0/ideal_AEP)
 
     cfg = nothing
@@ -86,13 +82,40 @@ function build_wind_farm_struct(x,turbine_x,turbine_y,turbine_z,hub_height,turbi
     opt_yaw && (turbine_yaw = Vector{input_type}(turbine_yaw))
     opt_diam && (rotor_diameter = Vector{input_type}(rotor_diameter))
 
-    preallocations = preallocations_struct(zeros(input_type,n_turbines,n_threads),zeros(input_type,n_turbines,n_threads),
-                    zeros(input_type,n_turbines,n_threads),zeros(input_type,n_turbines,n_threads),zeros(input_type,n_turbines,
-                    n_turbines,n_threads),zeros(input_type,n_turbines,n_turbines,n_threads),
-                    zeros(input_type,n_turbines,n_turbines,n_threads),zeros(input_type,n_turbines,n_turbines,n_threads))
+    preallocations = create_preallocations(turbine_x, turbine_y, turbine_yaw, rotor_diameter, hub_height, wind_farm_constants.wind_resource, wind_farm_constants.model_set)
 
     return wind_farm_struct(turbine_x, turbine_y, hub_height, turbine_yaw, rotor_diameter, results,
-                wind_farm_constants, AEP_scale, ideal_AEP, preallocations, update_function, AEP_gradient, AEP, cfg, force_single_thread)
+                wind_farm_constants, AEP_scale, ideal_AEP, preallocations, update_function, AEP_gradient, AEP, cfg)
+end
+
+function create_preallocations(turbine_x, turbine_y, turbine_yaw, rotor_diameter, hub_height, wind_resource, model_set)
+    n_turbines = length(turbine_x)
+    n_threads = Threads.nthreads()
+    n_states = determine_number_of_states(wind_resource, model_set)
+    T = promote_type(eltype(turbine_x),eltype(turbine_y),eltype(turbine_yaw),eltype(rotor_diameter),eltype(hub_height))
+    return preallocations_struct(zeros(T,n_turbines,n_threads),zeros(T,n_turbines,n_threads),
+                    zeros(T,n_turbines,n_threads),zeros(T,n_turbines,n_threads),zeros(T,n_turbines,
+                    n_turbines,n_threads),zeros(T,n_turbines,n_turbines,n_threads),
+                    zeros(T,n_turbines,n_turbines,n_threads),zeros(T,n_turbines,n_turbines,n_threads),
+                    zeros(T,n_turbines,n_threads),zeros(T,n_turbines,n_threads),
+                    zeros(T,n_turbines,n_threads),
+                    zeros(Int,n_turbines,n_threads),
+                    zeros(T,n_turbines,n_threads),
+                    zeros(T,n_turbines,n_threads),
+                    zeros(T,n_turbines,n_threads),
+                    zeros(T,n_states))
+end
+
+function determine_number_of_states(wind_resource, model_set)
+    if typeof(model_set.wake_combination_model) == SumOfSquaresFreestreamSuperposition
+        # find unique directions
+        unique_directions = unique(wind_resource.wind_directions)
+
+        # find how many unique directions there are
+        return length(unique_directions)
+    else
+        return length(wind_resource.wind_directions)
+    end
 end
 
 """
@@ -171,6 +194,21 @@ function calculate_aep!(farm,x)
 - `farm`: The wind_farm_struct
 - `x`: Vector containing the design variables
 """
+function calculate_aep!(farm,x::Array{ForwardDiff.Dual{T,V,N}}) where {T,V,N}
+    farm.update_function(farm,x)
+
+    AEP = calculate_aep(farm.turbine_x, farm.turbine_y, farm.constants.turbine_z, farm.rotor_diameter,
+                farm.hub_height, farm.turbine_yaw, farm.constants.ct_models, farm.constants.generator_efficiency,
+                farm.constants.cut_in_speed, farm.constants.cut_out_speed, farm.constants.rated_speed,
+                farm.constants.rated_power, farm.constants.wind_resource, farm.constants.power_models,
+                farm.constants.model_set,rotor_sample_points_y=farm.constants.rotor_sample_points_y,
+                rotor_sample_points_z=farm.constants.rotor_sample_points_z,
+                preallocations = farm.preallocations
+                ) .* farm.AEP_scale
+
+    return AEP
+end
+
 function calculate_aep!(farm,x)
     farm.update_function(farm,x)
 
@@ -180,17 +218,11 @@ function calculate_aep!(farm,x)
                 farm.constants.rated_power, farm.constants.wind_resource, farm.constants.power_models,
                 farm.constants.model_set,rotor_sample_points_y=farm.constants.rotor_sample_points_y,
                 rotor_sample_points_z=farm.constants.rotor_sample_points_z,
-                prealloc_turbine_velocities=farm.preallocations.prealloc_turbine_velocities,
-                prealloc_turbine_ct=farm.preallocations.prealloc_turbine_ct,
-                prealloc_turbine_ai=farm.preallocations.prealloc_turbine_ai,
-                prealloc_turbine_local_ti=farm.preallocations.prealloc_turbine_local_ti,
-                prealloc_wake_deficits=farm.preallocations.prealloc_wake_deficits,
-                prealloc_contribution_matrix=farm.preallocations.prealloc_contribution_matrix,
-                prealloc_deflections=farm.preallocations.prealloc_deflections,
-                prealloc_sigma_squared=farm.preallocations.prealloc_sigma_squared,
-                force_single_thread=farm.force_single_thread
+                preallocations = farm.preallocations
                 ) .* farm.AEP_scale
 
+    farm.AEP[1] = AEP
+
     return AEP
 end
 
@@ -211,20 +243,6 @@ function calculate_aep_gradient!(farm,x)
     return farm.AEP[1], farm.AEP_gradient
 end
 
-"""
-calculate_aep_gradient!(farm,x,sparse_struct::T)
-
-function to calculate the AEP and its gradient for the wind farm using sparse methods
-
-# Arguments
-- `farm`: The wind_farm_struct
-- `x`: Vector containing the design variables
-"""
-function calculate_aep_gradient!(farm,x,sparse_struct::T) where T <: AbstractSparseMethod
-    calculate_aep_gradient!(farm,x,sparse_struct)
-    return farm.AEP[1], farm.AEP_gradient
-end
-
 """
 calculate_spacing!(spacing_vec,x,spacing_struct)
 

src/power_models.jl

@@ -48,6 +48,31 @@ function rotate_to_wind_direction(xlocs::AbstractArray, ylocs::AbstractArray, wi
     return x_cart, y_cart
 end
 
+function rotate_to_wind_direction!(x_cart, y_cart, xlocs::AbstractArray, ylocs::AbstractArray, wind_direction_met::Number; center=0.0)
+    # use radians
+
+    # convert from meteorological polar system (CW, 0 rad.=N) to standard polar system (CCW, 0 rad.=E)
+    wind_direction_cart = met2cart(wind_direction_met)
+
+    cos_wdr = cos(-wind_direction_cart)
+    sin_wdr = sin(-wind_direction_cart)
+
+    if center != 0.0
+        center_x = center[1]
+        center_y = center[2]
+    else
+        center_x = 0.0
+        center_y = 0.0
+    end
+    # convert to cartesian coordinates with wind to positive x
+    for i in eachindex(xlocs)
+        x_cart[i] = ((xlocs[i] - center_x) * cos_wdr - (ylocs[i] - center_y) * sin_wdr) + center_x
+        y_cart[i] = ((xlocs[i] - center_x) * sin_wdr + (ylocs[i] - center_y) * cos_wdr) + center_y
+    end
+
+    return x_cart, y_cart
+end
+
 function rotate_to_wind_direction(xlocs::Number, ylocs::Number, wind_direction_met::Number; center=[0.0,0.0])
     # use radians
 
@@ -1683,3 +1708,33 @@ function star_boundary(n, ri, ro, rotation=0.0)
     # return
     return round.(vertices, digits=6)
 end
+
+# Nonallocating sortperm function
+function quicksort!(idx, x, lo::Int, hi::Int)
+    if lo < hi
+        p = partition!(idx, x, lo, hi)
+        quicksort!(idx, x, lo, p - 1)
+        quicksort!(idx, x, p + 1, hi)
+    end
+end
+
+function partition!(idx, x, lo::Int, hi::Int)
+    pivot = x[idx[hi]]
+    i = lo - 1
+    @inbounds for j in lo:hi-1
+        if x[idx[j]] ≤ pivot
+            i += 1
+            idx[i], idx[j] = idx[j], idx[i]
+        end
+    end
+    idx[i+1], idx[hi] = idx[hi], idx[i+1]
+    return i + 1
+end
+
+function sortperm!(idx::AbstractVector{Int}, x::AbstractVector)
+    @inbounds for i in eachindex(idx)
+        idx[i] = i
+    end
+    quicksort!(idx, x, 1, length(idx))
+    return nothing
+end

src/sparsity_functions.jl

@@ -571,12 +571,14 @@ end
 
 function _gauss_yaw_discontinuity(dt, x0, ky, kz, yaw, ct)
     # for clarity, break out the terms in the equation
-    a = ky + kz*cos(yaw)
-    b = 4.0 * ky * kz * cos(yaw)*(ct - 1.0)
-    c = 2.0 * sqrt(8.0) * ky * kz
+    cy = cos(yaw)
+    a = ky + kz*cy
+    kykz4 = 4.0 * ky * kz
+    b = kykz4 * cy*(ct - 1.0)
+    c = kykz4 * sqrt(2)
 
     # distance from rotor to the last point where the wake model is undefined
-    discontinuity_point = x0 + dt * (a - sqrt(a^2 - b))/c
+    discontinuity_point = x0 + dt * (a - sqrt(a*a - b))/c
 
     return discontinuity_point
 end
@@ -599,10 +601,20 @@ function _gauss_yaw_model_deficit(dx, dy, dz, dt, yaw, ct, ti, as, bs, ky, kz, w
         sigma_z = _gauss_yaw_spread_interpolated(dt, kz, dx, x0, 0.0, xd)
 
         # calculate velocity deficit #check - infty when large input ~= 500
-        ey = exp(-0.5*(dy/(wf*sigma_y))^2)
-        ez = exp(-0.5*(dz/(wf*sigma_z))^2)
+        c1 = -0.5*(dy/(wf*sigma_y))^2
+        c2 = -0.5*(dz/(wf*sigma_z))^2
 
-        sqrtterm = 1.0-ct*cos(yaw)/(8.0*(sigma_y*sigma_z/dt^2))
+        if c1 <= -750.0
+            loss = 0.0
+        end
+        if c2 <= -750.0
+            loss = 0.0
+        end
+        ey = exp(c1)
+        ez = exp(c2)
+
+        dt2 = dt*dt
+        sqrtterm = 1.0-ct*cos(yaw)/(8.0*(sigma_y*sigma_z/dt2))
         if sqrtterm >= 1e-8 #check - could try increasing tolerance
             loss = (1.0-sqrt(sqrtterm))*ey*ez
         else

src/user_functions.jl

@@ -182,10 +182,21 @@ function _bpa_theta_0(yaw, ct)
 end
 
 function _bpa_deflection(diam, ct, yaw, ky, kz, sigmay, sigmaz, theta0, x0)
-    a = theta0*x0/diam
-    b = (theta0/14.7)*sqrt(cos(yaw)/(ky*kz*ct))*(2.9-1.3*sqrt(1.0-ct)-ct)
-    c = (1.6+sqrt(ct))*(1.6*sqrt(8.0*sigmay*sigmaz/(cos(yaw)*diam^2))-ct)
-    d = (1.6-sqrt(ct))*(1.6*sqrt(8.0*sigmay*sigmaz/(cos(yaw)*diam^2))+ct)
+    diam_inv = 1.0/diam
+    a = theta0*x0*diam_inv
+    b1 = cos(yaw)
+    b = (theta0/14.7)*sqrt(b1/(ky*kz*ct))*(2.9-1.3*sqrt(1.0-ct)-ct)
+
+    c1 = 1.6*sqrt(8.0*sigmay*sigmaz * diam_inv * diam_inv / b1)
+    c2 = sqrt(ct)
+    
+    # c = (1.6+c2)*(c1-ct)
+    # d = (1.6-c2)*(c1+ct)
+    g1 = 1.6*c1-c2*ct
+    g2 = -1.6*ct+c1*c2
+    c = g1+g2
+    d = g1-g2
+
     y_deflection = diam*(a+b*log(c/d))
     return y_deflection
 end

src/utilities.jl

@@ -21,9 +21,8 @@ Unifying struct defining a wind farm and all necessary variables to calculate th
 - `AEP_gradient`: The gradient of the AEP
 - `AEP`: The AEP of the farm
 - `config`: The ForwardDiff config object if using ForwardDiff for AEP gradient calculation, otherwise nothing
-- `force_single_thread`: Boolean that forces the code to run in a single thread
 """
-struct wind_farm_struct{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15} <: AbstractWindFarmModel
+struct wind_farm_struct{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14} <: AbstractWindFarmModel
     turbine_x::T1
     turbine_y::T2
     hub_height::T3
@@ -38,7 +37,6 @@ struct wind_farm_struct{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15} <: A
     AEP_gradient::T12
     AEP::T13
     config::T14
-    force_single_thread::T15
 end
 
 """
@@ -56,7 +54,7 @@ struct that holds all the preallocated space for AEP calculation with one per th
 - `prealloc_deflections`: Matrix containing preallocated space for deflections
 - `prealloc_sigma_squared`: Matrix containing preallocated space for sigma squared
 """
-struct preallocations_struct{V,M}
+struct preallocations_struct{V,M,VI,Ve}
     prealloc_turbine_velocities::V
     prealloc_turbine_ct::V
     prealloc_turbine_ai::V
@@ -65,6 +63,14 @@ struct preallocations_struct{V,M}
     prealloc_contribution_matrix::M
     prealloc_deflections::M
     prealloc_sigma_squared::M
+    prealloc_rot_x::V
+    prealloc_rot_y::V
+    prealloc_power::V
+    prealloc_sort_index::VI
+    prealloc_diam::V
+    prealloc_hub::V
+    prealloc_yaw::V
+    prealloc_state_aep::Ve
 end
 
 """

src/wake_deficit_models.jl

@@ -22,7 +22,7 @@ function plot_velocity_row_deficit()
     include("./model_sets/model_set_2.jl")
 
     # calculate wind turbine velocities and corresponding aerodynamic operational states
-    turbine_inflow_velcities, turbine_ct, turbine_ai, turbine_local_ti = ff.turbine_velocities_one_direction(turbine_x, turbine_y, turbine_z, rotor_diameter, hub_height, turbine_yaw,
+    turbine_inflow_velcities, turbine_ct, turbine_ai, turbine_local_ti = ff.turbine_velocities_one_direction_full(turbine_x, turbine_y, turbine_z, rotor_diameter, hub_height, turbine_yaw,
     sorted_turbine_index, ct_model, rotor_sample_points_y, rotor_sample_points_z, windresource,
     model_set, wind_farm_state_id=1)