Remove interpolation allocation

examples/bench.jl

@@ -1,6 +1,8 @@
 using LinearAlgebra
 using ControlPlots
 using VortexStepMethod
+using Interpolations
+using BenchmarkTools
 
 plot = true
 
@@ -16,15 +18,40 @@ alpha = deg2rad(alpha_deg)
 # Step 2: Create wing geometry with linear panel distribution
 wing = Wing(n_panels, spanwise_panel_distribution=:linear)
 
-# Add wing sections - defining only tip sections with inviscid airfoil model
+# Test data generation
+n_angles = 5
+n_flaps = 5
+
+# Define angle ranges
+alphas = range(-deg2rad(10), deg2rad(10), n_angles)  # AoA range
+d_trailing_edge_angles = range(-deg2rad(30), deg2rad(30), n_flaps)  # Flap deflection range
+
+# Initialize coefficient matrices
+cl_matrix = zeros(n_angles, n_flaps)
+cd_matrix = zeros(n_angles, n_flaps)
+cm_matrix = zeros(n_angles, n_flaps)
+
+# Fill matrices with realistic aerodynamic data
+for (i, α) in enumerate(alphas)
+    for (j, δ) in enumerate(d_trailing_edge_angles)
+        cl_matrix[i,j] = 2π * α + π/2 * δ
+        cd_matrix[i,j] = 0.01 + 0.05 * α^2 + 0.03 * δ^2
+        cm_matrix[i,j] = -0.1 * α - 0.2 * δ
+    end
+end
+
+cl_interp = extrapolate(scale(interpolate(cl_matrix, BSpline(Linear())), alphas, d_trailing_edge_angles), NaN)
+cd_interp = extrapolate(scale(interpolate(cd_matrix, BSpline(Linear())), alphas, d_trailing_edge_angles), NaN)
+cm_interp = extrapolate(scale(interpolate(cm_matrix, BSpline(Linear())), alphas, d_trailing_edge_angles), NaN)
+
 add_section!(wing, 
     [0.0, span/2, 0.0],   # Left tip LE 
     [chord, span/2, 0.0],  # Left tip TE
-    :inviscid)
+    (:interpolations, (cl_interp, cd_interp, cm_interp)))
 add_section!(wing, 
     [0.0, -span/2, 0.0],  # Right tip LE
     [chord, -span/2, 0.0], # Right tip TE
-    :inviscid)
+    (:interpolations, (cl_interp, cd_interp, cm_interp)))
 
 # Step 3: Initialize aerodynamics
 wa = BodyAerodynamics([wing])
@@ -38,8 +65,9 @@ llt_solver = Solver(aerodynamic_model_type=:LLT)
 vsm_solver = Solver(aerodynamic_model_type=:VSM)
 
 # Step 5: Solve using both methods
-results_vsm = solve(vsm_solver, wa)
-@time results_vsm = solve(vsm_solver, wa)
+# @btime results_llt = solve($vsm_solver, $wa; log=false)
+# @btime results_vsm = solve($vsm_solver, $wa; log=false)
+@time results_vsm = solve(vsm_solver, wa; log=false)
 # time Python: 32.0 ms Ryzen 7950x
 # time Julia:   0.6 ms Ryzen 7950x
 #               0.8 ms laptop, performance mode, battery 

examples/ram_air_kite.jl

@@ -9,7 +9,7 @@ end
 using CSV
 using DataFrames
 
-plot = true
+plot = false
 
 # Create wing geometry
 wing = KiteWing("data/ram_air_kite_body.obj", "data/ram_air_kite_foil.dat")

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 = SMatrix{size(cl_matrix)...}(cl_matrix)
+    cd_matrix = SMatrix{size(cd_matrix)...}(cd_matrix)
+    cm_matrix = SMatrix{size(cm_matrix)...}(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::SMatrix, cd_matrix::SMatrix, cm_matrix::SMatrix) = 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,8 @@
 using LinearAlgebra
 
+const Interp1 = 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 Interp2 = Interpolations.FilledExtrapolation{Float64, 2, Interpolations.ScaledInterpolation{Float64, 2, Interpolations.BSplineInterpolation{Float64, 2, Interpolations.OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, Interpolations.BSpline{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Base.Slice{UnitRange{Int64}}, Base.Slice{UnitRange{Int64}}}}, Interpolations.BSpline{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, Interpolations.BSpline{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Float64}
+
 """
     Panel
 
@@ -36,9 +39,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, Interp1, Interp2}
+    cd_interp::Union{Nothing, Interp1, Interp2}
+    cm_interp::Union{Nothing, Interp1, Interp2}
     aerodynamic_center::MVec3
     control_point::MVec3
     bound_point_1::MVec3
@@ -82,45 +85,58 @@ 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)
+        ]
+
         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 = (
+                    (aero_1[2] + aero_2[2]) / 2,
+                    (aero_1[3] + aero_2[3]) / 2,
+                    (aero_1[4] + aero_2[4]) / 2
+                )
+                cl_interp = linear_interpolation(aero_1[1], polar_data[1]; extrapolation_bc=NaN)
+                cd_interp = linear_interpolation(aero_1[1], polar_data[2]; extrapolation_bc=NaN)
+                cm_interp = linear_interpolation(aero_1[1], 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 = (
+                    (aero_1[3] + aero_2[3]) / 2,
+                    (aero_1[4] + aero_2[4]) / 2,
+                    (aero_1[5] + aero_2[5]) / 2
+                )
+                cl_interp = linear_interpolation((aero_1[1], aero_1[2]), polar_data[1]; extrapolation_bc=NaN)
+                cd_interp = linear_interpolation((aero_1[1], aero_1[2]), polar_data[2]; extrapolation_bc=NaN)
+                cm_interp = linear_interpolation((aero_1[1], aero_1[2]), polar_data[3]; extrapolation_bc=NaN)
+                # cl_interp = extrapolate(scale(interpolate(cl_matrix, BSpline(Linear())), alphas, d_trailing_edge_angles), NaN)
+                # cd_interp = extrapolate(scale(interpolate(cl_matrix, BSpline(Linear())), alphas, d_trailing_edge_angles), NaN)
+                # cm_interp = extrapolate(scale(interpolate(cl_matrix, BSpline(Linear())), alphas, d_trailing_edge_angles), NaN)
+            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 +278,12 @@ 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, Interp1)
+            cl = panel.cl_interp(alpha)::Float64
+        elseif isa(panel.cl_interp, Interp2)
+            cl = panel.cl_interp(alpha, 0.0)::Float64
+        end
     else
         throw(ArgumentError("Unsupported aero model: $(panel.aero_model)"))
     end
@@ -284,9 +304,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, Interp1)
+            cd = panel.cd_interp(alpha)::Float64
+            cm = panel.cm_interp(alpha)::Float64
+        elseif isa(panel.cd_interp, Interp2)
+            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

@@ -127,7 +127,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,
@@ -138,7 +138,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