Remove interpolation allocation (#67)

* make example bench more proper

* working low allocation ram kite

* fix tests

* remove global

* add comment

* change back to old bench

* change center name

* remove second parameter syntax

Author: 1-Bart-1

examples/ram_air_kite.jl

@@ -20,10 +20,6 @@ VSM = Solver(
     aerodynamic_model_type=:VSM,
     is_with_artificial_damping=false
 )
-VSM_with_stall_correction = Solver(
-    aerodynamic_model_type=:VSM,
-    is_with_artificial_damping=true
-)
 
 # Setting velocity conditions
 v_a = 15.0

scripts/polars.jl

@@ -1,4 +1,4 @@
-using Distributed, Timers, Serialization, SharedArrays
+using Distributed, Timers, Serialization, SharedArrays, StaticArrays
 using Interpolations
 using Xfoil
 using ControlPlots
@@ -204,6 +204,10 @@ try
         cl_matrix[:, j], cd_matrix[:, j], cm_matrix[:, j] = run_solve_alpha(alphas, d_trailing_edge_angles[j], 
                         reynolds_number, x, y, lower, upper, kite_speed, SPEED_OF_SOUND, x_turn)
     end
+    cl_matrix = Matrix{Float64}(cl_matrix)
+    cd_matrix = Matrix{Float64}(cd_matrix)
+    cm_matrix = Matrix{Float64}(cm_matrix)
+
     display(cl_matrix)
 
     function plot_values(alphas, d_trailing_edge_angles, matrix, interp, name)
@@ -213,7 +217,8 @@ try
         X_data = collect(d_trailing_edge_angles) .+ zeros(length(alphas))'
         Y_data = collect(alphas)' .+ zeros(length(d_trailing_edge_angles))
 
-        interp_matrix = similar(matrix)
+        matrix = Matrix{Float64}(matrix)
+        interp_matrix = zeros(size(matrix)...)
         int_alphas, int_d_trailing_edge_angles = alphas .+ deg2rad(0.5), d_trailing_edge_angles .+ deg2rad(0.5)
         interp_matrix .= [interp(alpha, d_trailing_edge_angle) for alpha in int_alphas, d_trailing_edge_angle in int_d_trailing_edge_angles]
         X_int = collect(int_d_trailing_edge_angles) .+ zeros(length(int_alphas))'
@@ -241,7 +246,7 @@ try
     @info "Relative trailing_edge height: $(upper - lower)"
     @info "Reynolds number for flying speed of $kite_speed is $reynolds_number"
 
-    serialize(polar_path, (cl_interp, cd_interp, cm_interp))
+    serialize(polar_path, (alphas, d_trailing_edge_angles, cl_matrix, cd_matrix, cm_matrix))
 
 catch e
     @info "Removing processes"

src/VortexStepMethod.jl

@@ -10,7 +10,8 @@ using ControlPlots
 using Measures
 using LaTeXStrings
 using NonlinearSolve
-using Interpolations: linear_interpolation, Line, Extrapolation
+using Interpolations
+using Interpolations: linear_interpolation, Line, Extrapolation, FilledExtrapolation
 using Serialization
 using SharedArrays
 

src/kite_geometry.jl

@@ -241,15 +241,15 @@ mutable struct KiteWing <: AbstractWing
         end
 
         @info "Loding polars and kite info from $polar_path and $info_path"
-        (cl_interp, cd_interp, cm_interp) = deserialize(polar_path)
+        (alpha_range, beta_range, cl_matrix::Matrix, cd_matrix::Matrix, cm_matrix::Matrix) = deserialize(polar_path)
 
         (center_of_mass, inertia_tensor, circle_center_z, radius, gamma_tip, 
             le_interp, te_interp, area_interp) = deserialize(info_path)
 
         # Create sections
         sections = Section[]
         for gamma in range(-gamma_tip, gamma_tip, n_sections)
-            aero_input = (:interpolations, (cl_interp, cd_interp, cm_interp))
+            aero_input = (:polar_data, (alpha_range, beta_range, cl_matrix, cd_matrix, cm_matrix))
             LE_point = [0.0, 0.0, circle_center_z] .+ [le_interp(gamma), sin(gamma) * radius, cos(gamma) * radius]
             if !isapprox(alpha, 0.0)
                 local_y_vec = [0.0, sin(-gamma), cos(gamma)] × [1.0, 0.0, 0.0]

src/panel.jl

@@ -1,5 +1,9 @@
 using LinearAlgebra
 
+# static types for interpolations
+const I1 = Interpolations.FilledExtrapolation{Float64, 1, Interpolations.GriddedInterpolation{Float64, 1, Vector{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}}}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Float64}
+const I2 = Interpolations.FilledExtrapolation{Float64, 2, Interpolations.GriddedInterpolation{Float64, 2, Matrix{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Float64}
+
 """
     Panel
 
@@ -36,9 +40,9 @@ mutable struct Panel
     cl_coefficients::Vector{Float64}
     cd_coefficients::Vector{Float64}
     cm_coefficients::Vector{Float64}
-    cl_interp::Function
-    cd_interp::Function
-    cm_interp::Function
+    cl_interp::Union{Nothing, I1, I2}
+    cd_interp::Union{Nothing, I1, I2}
+    cm_interp::Union{Nothing, I1, I2}
     aero_center::MVec3
     control_point::MVec3
     bound_point_1::MVec3
@@ -82,45 +86,65 @@ mutable struct Panel
 
         # Initialize aerodynamic properties
         cl_coeffs, cd_coeffs, cm_coeffs = zeros(3), zeros(3), zeros(3)
-        cl_interp, cd_interp, cm_interp = (α::Float64)->0.0, (α::Float64)->0.0, (α::Float64)->0.0
 
+        # Calculate width
+        width = norm(bound_point_2 - bound_point_1)
+        
+        # Initialize filaments
+        filaments = [
+            BoundFilament(bound_point_2, bound_point_1),
+            BoundFilament(bound_point_1, TE_point_1),
+            BoundFilament(TE_point_2, bound_point_2)
+        ]
+
+        cl_interp, cd_interp, cm_interp = nothing, nothing, nothing
+
         if aero_model === :lei_airfoil_breukels
             cl_coeffs, cd_coeffs, cm_coeffs = compute_lei_coefficients(section_1, section_2)
+
         elseif aero_model === :polar_data
             aero_1 = section_1.aero_input[2]
             aero_2 = section_2.aero_input[2]
-            if size(aero_1) != size(aero_2)
+            if !all(size.(aero_1) .== size.(aero_2))
                 throw(ArgumentError("Polar data must have same shape"))
             end
-            polar_data = (aero_1 + aero_2) / 2
-            cl_interp_ = linear_interpolation(polar_data[:,1], polar_data[:,2]; extrapolation_bc=NaN)
-            cd_interp_ = linear_interpolation(polar_data[:,1], polar_data[:,3]; extrapolation_bc=NaN)
-            cm_interp_ = linear_interpolation(polar_data[:,1], polar_data[:,4]; extrapolation_bc=NaN)
-            cl_interp = (α::Float64) -> cl_interp_(α)
-            cd_interp = (α::Float64) -> cd_interp_(α)
-            cm_interp = (α::Float64) -> cm_interp_(α)
-        elseif aero_model === :interpolations
-            cl_left, cd_left, cm_left = section_1.aero_input[2]
-            cl_right, cd_right, cm_right = section_2.aero_input[2]
-            cl_interp = cl_left === cl_right ? (α::Float64) -> cl_left(α, 0.0) :
-                (α::Float64) -> 0.5cl_left(α, 0.0) + 0.5cl_right(α, 0.0) # TODO: add trailing edge deflection
-            cd_interp = cd_left === cd_right ? (α::Float64) -> cd_left(α, 0.0) :
-                (α::Float64) -> 0.5cd_left(α, 0.0) + 0.5cd_right(α, 0.0)
-            cm_interp = cm_left === cm_right ? (α::Float64) -> cm_left(α, 0.0) :
-                (α::Float64) -> 0.5cm_left(α, 0.0) + 0.5cm_right(α, 0.0)
+
+            if length(aero_1) == 4
+                !all(isapprox.(aero_1[1], aero_2[1])) && @error "Make sure you use the same alpha range for all your interpolations."
+
+                polar_data = (
+                    Vector{Float64}((aero_1[2] + aero_2[2]) / 2),
+                    Vector{Float64}((aero_1[3] + aero_2[3]) / 2),
+                    Vector{Float64}((aero_1[4] + aero_2[4]) / 2)
+                )
+                alphas = Vector{Float64}(aero_1[1])
+
+                cl_interp = linear_interpolation(alphas, polar_data[1]; extrapolation_bc=NaN)
+                cd_interp = linear_interpolation(alphas, polar_data[2]; extrapolation_bc=NaN)
+                cm_interp = linear_interpolation(alphas, polar_data[3]; extrapolation_bc=NaN)
+
+            elseif length(aero_1) == 5
+                !all(isapprox.(aero_1[1], aero_2[1])) && @error "Make sure you use the same alpha range for all your interpolations."
+                !all(isapprox.(aero_1[2], aero_2[2])) && @error "Make sure you use the same beta range for all your interpolations."
+
+                polar_data = (
+                    Matrix{Float64}((aero_1[3] + aero_2[3]) / 2),
+                    Matrix{Float64}((aero_1[4] + aero_2[4]) / 2),
+                    Matrix{Float64}((aero_1[5] + aero_2[5]) / 2)
+                )
+                alphas = Vector{Float64}(aero_1[1])
+                betas = Vector{Float64}(aero_1[2])
+
+                cl_interp = linear_interpolation((alphas, betas), polar_data[1]; extrapolation_bc=NaN)
+                cd_interp = linear_interpolation((alphas, betas), polar_data[2]; extrapolation_bc=NaN)
+                cm_interp = linear_interpolation((alphas, betas), polar_data[3]; extrapolation_bc=NaN)
+            else
+                throw(ArgumentError("Polar data in wrong format: $aero_1"))
+            end
+
         elseif !(aero_model === :inviscid)
             throw(ArgumentError("Unsupported aero model: $aero_model"))
         end
-        
-        # Calculate width
-        width = norm(bound_point_2 - bound_point_1)
-        
-        # Initialize filaments
-        filaments = [
-            BoundFilament(bound_point_2, bound_point_1),
-            BoundFilament(bound_point_1, TE_point_1),
-            BoundFilament(TE_point_2, bound_point_2)
-        ]
 
         new(
             TE_point_1, LE_point_1, TE_point_2, LE_point_2,
@@ -262,8 +286,14 @@ function calculate_cl(panel::Panel, alpha::Float64)::Float64
         end
     elseif panel.aero_model === :inviscid
         cl = 2π * alpha
-    elseif panel.aero_model === :polar_data || panel.aero_model === :interpolations
-        cl = panel.cl_interp(alpha)::Float64
+    elseif panel.aero_model === :polar_data
+        if isa(panel.cl_interp, I1)
+            cl = panel.cl_interp(alpha)::Float64
+        elseif isa(panel.cl_interp, I2)
+            cl = panel.cl_interp(alpha, 0.0)::Float64
+        else
+            throw(ArgumentError("cl_interp is $(panel.cl_interp)"))
+        end
     else
         throw(ArgumentError("Unsupported aero model: $(panel.aero_model)"))
     end
@@ -284,9 +314,14 @@ function calculate_cd_cm(panel::Panel, alpha::Float64)
         if abs(alpha) > (π/9)  # Outside ±20 degrees
             cd = 2 * sin(alpha)^3
         end
-    elseif panel.aero_model === :polar_data || panel.aero_model === :interpolations
-        cd = panel.cd_interp(alpha)::Float64
-        cm = panel.cm_interp(alpha)::Float64
+    elseif panel.aero_model === :polar_data
+        if isa(panel.cd_interp, I1)
+            cd = panel.cd_interp(alpha)::Float64
+            cm = panel.cm_interp(alpha)::Float64
+        elseif isa(panel.cd_interp, I2)
+            cd = panel.cd_interp(alpha, 0.0)::Float64
+            cm = panel.cm_interp(alpha, 0.0)::Float64
+        end
     elseif !(panel.aero_model === :inviscid)
         throw(ArgumentError("Unsupported aero model: $(panel.aero_model)"))
     end

src/solver.jl

@@ -147,7 +147,7 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n
     )
     # Try again with reduced relaxation factor if not converged
     if !converged && relaxation_factor > 1e-3
-        @warn "Running again with half the relaxation_factor = $(relaxation_factor/2)"
+        log && @warn "Running again with half the relaxation_factor = $(relaxation_factor/2)"
         converged, gamma_new, alpha_array, v_a_array = gamma_loop(
             solver,
             body_aero,
@@ -158,7 +158,8 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n
             y_airf_array,
             z_airf_array,
             panels,
-            relaxation_factor/2
+            relaxation_factor/2;
+            log
         )
     end
 

src/wing_geometry.jl

@@ -68,7 +68,8 @@ Add a new section to the wing.
 - aero_input: Can be:
   - :inviscid
   - :`lei_airfoil_breukels`
-  - (:interpolations, (`cl_interp`, `cd_interp`, `cm_interp`))
+  - (:polar_data, (`alpha_range`, `cl_vector`, `cd_vector`, `cm_vector`))
+  - (:polar_data, (`alpha_range`, `beta_range`, `cl_matrix`, `cd_matrix`, `cm_matrix`))
 """
 function add_section!(wing::Wing, LE_point::Vector{Float64}, 
                      TE_point::Vector{Float64}, aero_input)
@@ -157,34 +158,6 @@ function refine_aerodynamic_mesh(wing::AbstractWing)
     end
 end
 
-"""
-    interpolate_to_common_alpha(alpha_common::Vector{Float64}, 
-                              alpha_orig::Vector{Float64},
-                              CL_orig::Vector{Float64},
-                              CD_orig::Vector{Float64},
-                              CM_orig::Vector{Float64})
-
-Interpolate aerodynamic coefficients to a common alpha range.
-"""
-function interpolate_to_common_alpha(alpha_common, 
-                                   alpha_orig,
-                                   CL_orig,
-                                   CD_orig,
-                                   CM_orig)
-
-    # Create interpolation objects
-    itp_CL = linear_interpolation(alpha_orig, CL_orig)
-    itp_CD = linear_interpolation(alpha_orig, CD_orig)
-    itp_CM = linear_interpolation(alpha_orig, CM_orig)
-    
-    # Evaluate at common alpha points
-    CL_common = itp_CL.(alpha_common)
-    CD_common = itp_CD.(alpha_common)
-    CM_common = itp_CM.(alpha_common)
-
-    return CL_common, CD_common, CM_common
-end
-
 """
     calculate_new_aero_input(aero_input, 
                             section_index::Int,
@@ -212,46 +185,38 @@ function calculate_new_aero_input(aero_input,
         polar_right = aero_input[section_index + 1][2]
 
         # Unpack polar data
-        @views begin
-            alpha_left, CL_left, CD_left, CM_left = (
-                polar_left[:, i] for i in 1:4
-            )
-            alpha_right, CL_right, CD_right, CM_right = (
-                polar_right[:, i] for i in 1:4
-            )
-        end
-        
-        # Create common alpha array
-        alpha_common = sort(unique(vcat(alpha_left, alpha_right)))
-        
-        # Interpolate both polars
-        CL_left_common, CD_left_common, CM_left_common = interpolate_to_common_alpha(
-            alpha_common, alpha_left, CL_left, CD_left, CM_left
-        )
-        CL_right_common, CD_right_common, CM_right_common = interpolate_to_common_alpha(
-            alpha_common, alpha_right, CL_right, CD_right, CM_right
-        )
-        
-        # Weighted interpolation
-        CL_interp = CL_left_common .* left_weight .+ CL_right_common .* right_weight
-        CD_interp = CD_left_common .* left_weight .+ CD_right_common .* right_weight
-        CM_interp = CM_left_common .* left_weight .+ CM_right_common .* right_weight
-        
-        return (:polar_data, hcat(alpha_common, CD_interp, CL_interp, CM_interp))
+        if length(polar_left) == 4
+            alpha_left, CL_left, CD_left, CM_left = polar_left
+            alpha_right, CL_right, CD_right, CM_right = polar_right
 
-    elseif model_type === :interpolations
-        cl_left, cd_left, cm_left = aero_input[section_index][2]
-        cl_right, cd_right, cm_right = aero_input[section_index + 1][2]
-
-        cl_interp = cl_left === cl_right ? (α, β) -> cl_left[α, β] :
-            (α, β) -> left_weight * cl_left[α, β] + right_weight * cl_right[α, β]
-        cd_interp = cd_left === cd_right ? (α, β) -> cd_left[α, β] :
-            (α, β) -> left_weight * cd_left[α, β] + right_weight * cd_right[α, β]
-        cm_interp = cm_left === cm_right ? (α, β) -> cm_left[α, β] :
-            (α, β) -> left_weight * cm_left[α, β] + right_weight * cm_right[α, β]
+            # Create common alpha array
+            !all(isapprox.(alpha_left, alpha_right)) && @error "Make sure you use the same alpha range for all your interpolations."
+            !isa(CL_right, AbstractVector) && @error "Provide polar data in the correct format: (alpha, cl, cd, cm)"
+            
+            # Weighted interpolation
+            CL_data = CL_left .* left_weight .+ CL_right .* right_weight
+            CD_data = CD_left .* left_weight .+ CD_right .* right_weight
+            CM_data = CM_left .* left_weight .+ CM_right .* right_weight
+            
+            return (:polar_data, (alpha_left, CL_data, CD_data, CM_data))
+            
+        elseif length(polar_left) == 5
+            alpha_left, beta_left, CL_left, CD_left, CM_left = polar_left
+            alpha_right, beta_right, CL_right, CD_right, CM_right = polar_right
+            
+            # Create common alpha array
+            !all(isapprox.(alpha_left, alpha_right)) && @error "Make sure you use the same alpha range for all your interpolations."
+            !all(isapprox.(beta_left, beta_right)) && @error "Make sure you use the same alpha range for all your interpolations."
+            !isa(CL_right, AbstractMatrix) && @error "Provide polar data in the correct format: (alpha, beta, cl, cd, cm)"
+
+            # Weighted interpolation
+            CL_data = CL_left .* left_weight .+ CL_right .* right_weight
+            CD_data = CD_left .* left_weight .+ CD_right .* right_weight
+            CM_data = CM_left .* left_weight .+ CM_right .* right_weight
+            
+            return (:polar_data, (alpha_left, beta_left, CL_data, CD_data, CM_data))
+        end
 
-        return (:interpolations, (cl_interp, cd_interp, cm_interp))
-        
     elseif model_type === :lei_airfoil_breukels
         tube_diameter_left = aero_input[section_index][2][1]
         tube_diameter_right = aero_input[section_index + 1][2][1]