From e72739d4f96ca9fe224708e56d11a6191942105c Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Mon, 24 Feb 2025 20:36:21 -0600 Subject: [PATCH 1/8] make example bench more proper --- examples/bench.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/bench.jl b/examples/bench.jl index 19e45634..0ef6757c 100644 --- a/examples/bench.jl +++ b/examples/bench.jl @@ -1,6 +1,7 @@ using LinearAlgebra using ControlPlots using VortexStepMethod +using BenchmarkTools plot = true @@ -38,8 +39,8 @@ 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) +@btime results_vsm = solve($vsm_solver, $wa) # time Python: 32.0 ms Ryzen 7950x # time Julia: 0.6 ms Ryzen 7950x # 0.8 ms laptop, performance mode, battery From 773ecd6044549a99a9e1664897eb641bb4ca863e Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Tue, 25 Feb 2025 13:40:17 -0600 Subject: [PATCH 2/8] working low allocation ram kite --- examples/bench.jl | 37 ++++++++++++--- examples/ram_air_kite.jl | 2 +- scripts/polars.jl | 11 +++-- src/VortexStepMethod.jl | 3 +- src/kite_geometry.jl | 4 +- src/panel.jl | 97 +++++++++++++++++++++++++--------------- src/solver.jl | 5 ++- src/wing_geometry.jl | 94 ++++++++++++-------------------------- test/bench.jl | 82 ++++++++++++++++++++------------- 9 files changed, 189 insertions(+), 146 deletions(-) diff --git a/examples/bench.jl b/examples/bench.jl index 0ef6757c..e39ae495 100644 --- a/examples/bench.jl +++ b/examples/bench.jl @@ -1,6 +1,7 @@ using LinearAlgebra using ControlPlots using VortexStepMethod +using Interpolations using BenchmarkTools plot = true @@ -17,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]) @@ -39,8 +65,9 @@ llt_solver = Solver(aerodynamic_model_type=:LLT) vsm_solver = Solver(aerodynamic_model_type=:VSM) # Step 5: Solve using both methods -@btime results_llt = solve($vsm_solver, $wa) -@btime 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 diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl index 3510a3df..2990e307 100644 --- a/examples/ram_air_kite.jl +++ b/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") diff --git a/scripts/polars.jl b/scripts/polars.jl index bf7c5d1e..0fb4458c 100644 --- a/scripts/polars.jl +++ b/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" diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 27483333..48e77cee 100644 --- a/src/VortexStepMethod.jl +++ b/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 diff --git a/src/kite_geometry.jl b/src/kite_geometry.jl index 51cb6d9f..6193d729 100644 --- a/src/kite_geometry.jl +++ b/src/kite_geometry.jl @@ -241,7 +241,7 @@ 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) @@ -249,7 +249,7 @@ mutable struct KiteWing <: AbstractWing # 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] diff --git a/src/panel.jl b/src/panel.jl index 9841c7b9..eeb32494 100644 --- a/src/panel.jl +++ b/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 diff --git a/src/solver.jl b/src/solver.jl index d2652763..9949e07b 100644 --- a/src/solver.jl +++ b/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 diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index fa21d422..aba6be50 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -157,34 +157,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 +184,36 @@ 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." + + # 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, hcat(alpha_left, CD_data, CL_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." + + # 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] diff --git a/test/bench.jl b/test/bench.jl index 03cbe84c..c2feda9b 100644 --- a/test/bench.jl +++ b/test/bench.jl @@ -33,7 +33,7 @@ using LinearAlgebra [chord, -span/2, 0.0], # Right tip TE :inviscid) - body_aero = BodyAerodynamics([wing]) + global body_aero = BodyAerodynamics([wing]) vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a set_va!(body_aero, (vel_app, 0.0)) # Second parameter is yaw rate @@ -60,7 +60,7 @@ using LinearAlgebra for model in models for frac in core_radius_fractions @testset "Model $model Core Radius Fraction $frac" begin - global result = @benchmark calculate_AIC_matrices!($body_aero, $model, $frac, $va_norm_array, $va_unit_array) samples = 100 + result = @benchmark calculate_AIC_matrices!($body_aero, $model, $frac, $va_norm_array, $va_unit_array) samples = 1 evals = 1 @test result.allocs ≤ 100 @info "Model: $(model) \t Core radius fraction: $(frac) \t Allocations: $(result.allocs) \t Memory: $(result.memory)" end @@ -86,25 +86,45 @@ using LinearAlgebra y_airf_array[i, :] .= panel.y_airf z_airf_array[i, :] .= panel.z_airf end + + n_angles = 5 + alphas = range(-deg2rad(10), deg2rad(10), n_angles) + cls = [2π * α for α in alphas] + cds = [0.01 + 0.05 * α^2 for α in alphas] + cms = [-0.1 * α for α in alphas] + for model in models - solver = Solver( - aerodynamic_model_type=model - ) - result = @benchmark gamma_loop( - $solver, - $body_aero, - $gamma_new, - $va_array, - $chord_array, - $x_airf_array, - $y_airf_array, - $z_airf_array, - $body_aero.panels, - 0.5; - log = false - ) samples = 100 - @test result.allocs ≤ 100 - @info "Model: $model \t Allocations: $(result.allocs) Memory: $(result.memory)" + for aero_model in [:inviscid, (:polar_data, [alphas, cls, cds, cms])] + wing = Wing(n_panels, spanwise_panel_distribution=:linear) + add_section!(wing, + [0.0, span/2, 0.0], # Left tip LE + [chord, span/2, 0.0], # Left tip TE + aero_model) + add_section!(wing, + [0.0, -span/2, 0.0], # Right tip LE + [chord, -span/2, 0.0], # Right tip TE + aero_model) + body_aero = BodyAerodynamics([wing]) + + solver = Solver( + aerodynamic_model_type=model + ) + result = @benchmark gamma_loop( + $solver, + $body_aero, + $gamma_new, + $va_array, + $chord_array, + $x_airf_array, + $y_airf_array, + $z_airf_array, + $body_aero.panels, + 0.5; + log = false + ) samples = 1 evals = 1 + @test result.allocs ≤ 100 + @info "Model: $model \t Aero_model: $aero_model \t Allocations: $(result.allocs) Memory: $(result.memory)" + end end end @@ -120,16 +140,18 @@ using LinearAlgebra va_norm_array = zeros(n_panels) va_unit_array = zeros(n_panels, 3) - # Fill arrays with data - for (i, panel) in enumerate(body_aero.panels) - chord_array[i] = panel.chord - x_airf_array[i, :] .= panel.x_airf - y_airf_array[i, :] .= panel.y_airf - z_airf_array[i, :] .= panel.z_airf - va_array[i, :] .= panel.va - end - + + # # Fill arrays with data + # for (i, panel) in enumerate(body_aero.panels) + # chord_array[i] = panel.chord + # x_airf_array[i, :] .= panel.x_airf + # y_airf_array[i, :] .= panel.y_airf + # z_airf_array[i, :] .= panel.z_airf + # va_array[i, :] .= panel.va + # end + set_va!(body_aero, (vel_app, 0.0)) # Second parameter is yaw rate results = @MVector zeros(3) + result = @benchmark calculate_results( $body_aero, $gamma, @@ -148,7 +170,7 @@ using LinearAlgebra $va_unit_array, $body_aero.panels, false - ) + ) samples = 1 evals = 1 @test_broken result.allocs ≤ 100 end From 9451b66de36ba153d0bf11e5fdb9efcf537a5870 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Tue, 25 Feb 2025 20:38:34 -0600 Subject: [PATCH 3/8] fix tests --- examples/bench.jl | 103 +++++++++++++------------------------ examples/ram_air_kite.jl | 6 +-- scripts/polars.jl | 6 +-- src/kite_geometry.jl | 2 +- src/panel.jl | 57 +++++++++++--------- src/wing_geometry.jl | 7 ++- test/bench.jl | 2 +- test/test_kite_geometry.jl | 7 +-- test/test_panel.jl | 20 ++++--- 9 files changed, 93 insertions(+), 117 deletions(-) diff --git a/examples/bench.jl b/examples/bench.jl index e39ae495..fa93aa93 100644 --- a/examples/bench.jl +++ b/examples/bench.jl @@ -1,75 +1,42 @@ -using LinearAlgebra -using ControlPlots using VortexStepMethod -using Interpolations -using BenchmarkTools - -plot = true - -# Step 1: Define wing parameters -n_panels = 20 # Number of panels -span = 20.0 # Wing span [m] -chord = 1.0 # Chord length [m] -v_a = 20.0 # Magnitude of inflow velocity [m/s] -density = 1.225 # Air density [kg/m³] -alpha_deg = 30.0 # Angle of attack [degrees] -alpha = deg2rad(alpha_deg) - -# Step 2: Create wing geometry with linear panel distribution -wing = Wing(n_panels, spanwise_panel_distribution=:linear) - -# 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) +using LinearAlgebra +using Pkg -# 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 +if !("CSV" ∈ keys(Pkg.project().dependencies)) + using TestEnv + TestEnv.activate() end +using CSV +using DataFrames +using BenchmarkTools -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 - (: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 - (:interpolations, (cl_interp, cd_interp, cm_interp))) - -# Step 3: Initialize aerodynamics -wa = BodyAerodynamics([wing]) - -# Set inflow conditions -vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a -set_va!(wa, (vel_app, 0.0)) # Second parameter is yaw rate - -# Step 4: Initialize solvers for both LLT and VSM methods -llt_solver = Solver(aerodynamic_model_type=:LLT) -vsm_solver = Solver(aerodynamic_model_type=:VSM) +plot = true -# Step 5: Solve using both methods -# @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 +# Create wing geometry +wing = KiteWing("data/ram_air_kite_body.obj", "data/ram_air_kite_foil.dat") +body_aero = BodyAerodynamics([wing]) + +# Create solvers +VSM = Solver( + aerodynamic_model_type=:VSM, + is_with_artificial_damping=false +) + +# Setting velocity conditions +v_a = 15.0 +aoa = 15.0 +side_slip = 0.0 +yaw_rate = 0.0 +aoa_rad = deg2rad(aoa) +vel_app = [ + cos(aoa_rad) * cos(side_slip), + sin(side_slip), + sin(aoa_rad) +] * v_a +body_aero.va = vel_app + +# Solving and plotting distributions +results = solve(VSM, body_aero) +@btime results = solve($VSM, $body_aero) nothing \ No newline at end of file diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl index 2990e307..250ea48b 100644 --- a/examples/ram_air_kite.jl +++ b/examples/ram_air_kite.jl @@ -9,7 +9,7 @@ end using CSV using DataFrames -plot = false +plot = true # Create wing geometry wing = KiteWing("data/ram_air_kite_body.obj", "data/ram_air_kite_foil.dat") @@ -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 diff --git a/scripts/polars.jl b/scripts/polars.jl index 0fb4458c..7f9de6ca 100644 --- a/scripts/polars.jl +++ b/scripts/polars.jl @@ -204,9 +204,9 @@ 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) + cl_matrix = Matrix{Float64}(cl_matrix) + cd_matrix = Matrix{Float64}(cd_matrix) + cm_matrix = Matrix{Float64}(cm_matrix) display(cl_matrix) diff --git a/src/kite_geometry.jl b/src/kite_geometry.jl index 6193d729..5c5c0665 100644 --- a/src/kite_geometry.jl +++ b/src/kite_geometry.jl @@ -241,7 +241,7 @@ mutable struct KiteWing <: AbstractWing end @info "Loding polars and kite info from $polar_path and $info_path" - (alpha_range, beta_range, cl_matrix::SMatrix, cd_matrix::SMatrix, cm_matrix::SMatrix) = 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) diff --git a/src/panel.jl b/src/panel.jl index eeb32494..28d09df8 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -1,7 +1,7 @@ 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} +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 @@ -39,9 +39,9 @@ mutable struct Panel cl_coefficients::Vector{Float64} cd_coefficients::Vector{Float64} cm_coefficients::Vector{Float64} - cl_interp::Union{Nothing, Interp1, Interp2} - cd_interp::Union{Nothing, Interp1, Interp2} - cm_interp::Union{Nothing, Interp1, Interp2} + cl_interp::Union{Nothing, I1, I2} + cd_interp::Union{Nothing, I1, I2} + cm_interp::Union{Nothing, I1, I2} aerodynamic_center::MVec3 control_point::MVec3 bound_point_1::MVec3 @@ -96,6 +96,8 @@ mutable struct Panel 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) @@ -108,30 +110,35 @@ mutable struct Panel 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 + 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) ) - 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) + 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 = ( - (aero_1[3] + aero_2[3]) / 2, - (aero_1[4] + aero_2[4]) / 2, - (aero_1[5] + aero_2[5]) / 2 + 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) ) - 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) + 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) @@ -279,10 +286,12 @@ function calculate_cl(panel::Panel, alpha::Float64)::Float64 elseif panel.aero_model === :inviscid cl = 2π * alpha elseif panel.aero_model === :polar_data - if isa(panel.cl_interp, Interp1) + if isa(panel.cl_interp, I1) cl = panel.cl_interp(alpha)::Float64 - elseif isa(panel.cl_interp, Interp2) + 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)")) @@ -305,10 +314,10 @@ function calculate_cd_cm(panel::Panel, alpha::Float64) cd = 2 * sin(alpha)^3 end elseif panel.aero_model === :polar_data - if isa(panel.cd_interp, Interp1) + if isa(panel.cd_interp, I1) cd = panel.cd_interp(alpha)::Float64 cm = panel.cm_interp(alpha)::Float64 - elseif isa(panel.cd_interp, Interp2) + elseif isa(panel.cd_interp, I2) cd = panel.cd_interp(alpha, 0.0)::Float64 cm = panel.cm_interp(alpha, 0.0)::Float64 end diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index aba6be50..c185897a 100644 --- a/src/wing_geometry.jl +++ b/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) @@ -190,13 +191,14 @@ function calculate_new_aero_input(aero_input, # 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, hcat(alpha_left, CD_data, CL_data, CM_data)) + 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 @@ -205,6 +207,7 @@ function calculate_new_aero_input(aero_input, # 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 diff --git a/test/bench.jl b/test/bench.jl index c2feda9b..77ecaca8 100644 --- a/test/bench.jl +++ b/test/bench.jl @@ -94,7 +94,7 @@ using LinearAlgebra cms = [-0.1 * α for α in alphas] for model in models - for aero_model in [:inviscid, (:polar_data, [alphas, cls, cds, cms])] + for aero_model in [:inviscid, (:polar_data, (alphas, cls, cds, cms))] wing = Wing(n_panels, spanwise_panel_distribution=:linear) add_section!(wing, [0.0, span/2, 0.0], # Left tip LE diff --git a/test/test_kite_geometry.jl b/test/test_kite_geometry.jl index bf962c2c..3b131853 100644 --- a/test/test_kite_geometry.jl +++ b/test/test_kite_geometry.jl @@ -100,14 +100,9 @@ using Serialization end end - # Create interpolations - 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) - # Serialize polar data polar_path = test_dat_path[1:end-4] * "_polar.bin" - serialize(polar_path, (cl_interp, cd_interp, cm_interp)) + serialize(polar_path, (alphas, d_trailing_edge_angles, cl_matrix, cd_matrix, cm_matrix)) # Create and serialize obj file faces = [[i, i+1, i+2] for i in 1:2:length(vertices)-2] diff --git a/test/test_panel.jl b/test/test_panel.jl index 7cbc795c..4b654488 100644 --- a/test/test_panel.jl +++ b/test/test_panel.jl @@ -115,16 +115,22 @@ end # Generate mock polar data alphas = deg2rad.(-10:1:25) n_points = length(alphas) - polar_data = zeros(n_points, 4) - big_polar_data = zeros(n_points, 4) + polar_data = [zeros(n_points) for _ in 1:4] + big_polar_data = [zeros(n_points) for _ in 1:4] # Fill polar data with realistic values for (i, alpha) in enumerate(alphas) cd = 0.015 + 0.015 * abs(alpha/10)^1.5 # Drag increases with angle cl = 0.1 * alpha + 0.01 * alpha^2 * exp(-alpha/20) # Lift with stall behavior cm = -0.02 * alpha # Linear pitching moment - polar_data[i,:] = [alpha, cl, cd, cm] - big_polar_data[i,:] = [alpha, 1.1cl, 1.1cd, 1.1cm] + polar_data[1][i] = alpha + polar_data[2][i] = cl + polar_data[3][i] = cd + polar_data[4][i] = cm + big_polar_data[1][i] = alpha + big_polar_data[2][i] = 1.1cl + big_polar_data[3][i] = 1.1cd + big_polar_data[4][i] = 1.1cm end # Create two sections with slightly different polar data @@ -142,9 +148,9 @@ end cd, cm = calculate_cd_cm(panel, alpha) # Calculate expected values through interpolation - expected_cl = linear_interpolation(polar_data[:,1], polar_data[:,2], extrapolation_bc=Line())(alpha) - expected_cd = linear_interpolation(polar_data[:,1], polar_data[:,3], extrapolation_bc=Line())(alpha) - expected_cm = linear_interpolation(polar_data[:,1], polar_data[:,4], extrapolation_bc=Line())(alpha) + expected_cl = linear_interpolation(polar_data[1], polar_data[2], extrapolation_bc=Line())(alpha) + expected_cd = linear_interpolation(polar_data[1], polar_data[3], extrapolation_bc=Line())(alpha) + expected_cm = linear_interpolation(polar_data[1], polar_data[4], extrapolation_bc=Line())(alpha) # Average with slightly different data (1.1 factor) expected_cl = (expected_cl + 1.1 * expected_cl) / 2 From 8aa1ca8c624d716c9880be64c3115fb780f77ddf Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Tue, 25 Feb 2025 20:40:11 -0600 Subject: [PATCH 4/8] remove global --- test/bench.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/bench.jl b/test/bench.jl index 77ecaca8..ddf0e6bb 100644 --- a/test/bench.jl +++ b/test/bench.jl @@ -33,7 +33,7 @@ using LinearAlgebra [chord, -span/2, 0.0], # Right tip TE :inviscid) - global body_aero = BodyAerodynamics([wing]) + body_aero = BodyAerodynamics([wing]) vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a set_va!(body_aero, (vel_app, 0.0)) # Second parameter is yaw rate From 554a4cf219c7f02f64dc2e4c40dadf239f5669dc Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Wed, 26 Feb 2025 08:24:00 -0600 Subject: [PATCH 5/8] add comment --- src/panel.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/panel.jl b/src/panel.jl index 28d09df8..c811d66b 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -1,5 +1,6 @@ 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} From 2c275c75a16bd003d6e4b2bc64056ec1779fc08e Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Wed, 26 Feb 2025 08:27:02 -0600 Subject: [PATCH 6/8] change back to old bench --- examples/bench.jl | 77 +++++++++++++++++++++++++---------------------- 1 file changed, 41 insertions(+), 36 deletions(-) diff --git a/examples/bench.jl b/examples/bench.jl index fa93aa93..19e45634 100644 --- a/examples/bench.jl +++ b/examples/bench.jl @@ -1,42 +1,47 @@ -using VortexStepMethod using LinearAlgebra -using Pkg - -if !("CSV" ∈ keys(Pkg.project().dependencies)) - using TestEnv - TestEnv.activate() -end -using CSV -using DataFrames -using BenchmarkTools +using ControlPlots +using VortexStepMethod plot = true -# Create wing geometry -wing = KiteWing("data/ram_air_kite_body.obj", "data/ram_air_kite_foil.dat") -body_aero = BodyAerodynamics([wing]) - -# Create solvers -VSM = Solver( - aerodynamic_model_type=:VSM, - is_with_artificial_damping=false -) - -# Setting velocity conditions -v_a = 15.0 -aoa = 15.0 -side_slip = 0.0 -yaw_rate = 0.0 -aoa_rad = deg2rad(aoa) -vel_app = [ - cos(aoa_rad) * cos(side_slip), - sin(side_slip), - sin(aoa_rad) -] * v_a -body_aero.va = vel_app - -# Solving and plotting distributions -results = solve(VSM, body_aero) -@btime results = solve($VSM, $body_aero) +# Step 1: Define wing parameters +n_panels = 20 # Number of panels +span = 20.0 # Wing span [m] +chord = 1.0 # Chord length [m] +v_a = 20.0 # Magnitude of inflow velocity [m/s] +density = 1.225 # Air density [kg/m³] +alpha_deg = 30.0 # Angle of attack [degrees] +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 +add_section!(wing, + [0.0, span/2, 0.0], # Left tip LE + [chord, span/2, 0.0], # Left tip TE + :inviscid) +add_section!(wing, + [0.0, -span/2, 0.0], # Right tip LE + [chord, -span/2, 0.0], # Right tip TE + :inviscid) + +# Step 3: Initialize aerodynamics +wa = BodyAerodynamics([wing]) + +# Set inflow conditions +vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a +set_va!(wa, (vel_app, 0.0)) # Second parameter is yaw rate + +# Step 4: Initialize solvers for both LLT and VSM methods +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) +# time Python: 32.0 ms Ryzen 7950x +# time Julia: 0.6 ms Ryzen 7950x +# 0.8 ms laptop, performance mode, battery nothing \ No newline at end of file From 7a8750d91434a6f16f8661482d1819004fc0b04d Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Wed, 26 Feb 2025 08:47:05 -0600 Subject: [PATCH 7/8] change center name --- src/panel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/panel.jl b/src/panel.jl index 29d30d9f..87db1ec9 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -43,7 +43,7 @@ mutable struct Panel cl_interp::Union{Nothing, I1, I2} cd_interp::Union{Nothing, I1, I2} cm_interp::Union{Nothing, I1, I2} - aerodynamic_center::MVec3 + aero_center::MVec3 control_point::MVec3 bound_point_1::MVec3 bound_point_2::MVec3 From aaad5a9fbab28be9e621c9d46985d796f4c166cf Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Wed, 26 Feb 2025 08:48:55 -0600 Subject: [PATCH 8/8] remove second parameter syntax --- test/bench.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/bench.jl b/test/bench.jl index dc313817..fdf3a459 100644 --- a/test/bench.jl +++ b/test/bench.jl @@ -150,7 +150,7 @@ using LinearAlgebra # z_airf_array[i, :] .= panel.z_airf # va_array[i, :] .= panel.va # end - set_va!(body_aero, (vel_app, 0.0)) # Second parameter is yaw rate + set_va!(body_aero, vel_app) results = @MVector zeros(3) result = @benchmark calculate_results(