From 9e2834a2ba6f741811ddb10d79d2a3ebf5050fff Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Thu, 20 Feb 2025 23:04:42 -0600 Subject: [PATCH 01/24] reduce allocations more --- src/solver.jl | 72 +++++++++++++++++++++++++++++++++++---------------- 1 file changed, 50 insertions(+), 22 deletions(-) diff --git a/src/solver.jl b/src/solver.jl index be540907..2723ce6c 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -74,7 +74,7 @@ function solve(solver::Solver, wing_aero::WingAerodynamics, gamma_distribution=n z_airf_array = zeros(n_panels, 3) va_array = zeros(n_panels, 3) chord_array = zeros(n_panels) - + # Fill arrays from panels for (i, panel) in enumerate(panels) x_airf_array[i, :] .= panel.x_airf @@ -170,6 +170,8 @@ function solve(solver::Solver, wing_aero::WingAerodynamics, gamma_distribution=n return results end +cross3(x,y) = cross(SVector{3,eltype(x)}(x), SVector{3,eltype(y)}(y)) + """ gamma_loop(solver::Solver, gamma_new::Vector{Float64}, AIC_x::Matrix{Float64}, AIC_y::Matrix{Float64}, AIC_z::Matrix{Float64}, va_array::Matrix{Float64}, @@ -195,33 +197,58 @@ function gamma_loop( relaxation_factor::Float64 ) converged = false + n_panels = wing_aero.n_panels alpha_array = wing_aero.alpha_array Umag_array = wing_aero.Umag_array - + Umagw_array = similar(Umag_array) + + gamma = copy(gamma_new) + abs_gamma_new = copy(gamma_new) + induced_velocity_all = zeros(size(AIC_x, 1), 3) + relative_velocity_array = similar(va_array) + relative_velocity_crossz = similar(relative_velocity_array) + Uinfcrossz_array = similar(va_array) + cl_array = zeros(n_panels) + damp = zeros(length(gamma)) + v_normal_array = zeros(n_panels) + v_tangential_array = zeros(n_panels) + iters = 0 for i in 1:solver.max_iterations iters += 1 - gamma = copy(gamma_new) + gamma .= gamma_new # Calculate induced velocities - induced_velocity_all = hcat( - AIC_x * gamma, - AIC_y * gamma, - AIC_z * gamma - ) + mul!(view(induced_velocity_all, :, 1), AIC_x, gamma) + mul!(view(induced_velocity_all, :, 2), AIC_y, gamma) + mul!(view(induced_velocity_all, :, 3), AIC_z, gamma) - relative_velocity_array = va_array .+ induced_velocity_all - relative_velocity_crossz = cross.(eachrow(relative_velocity_array), eachrow(z_airf_array)) - Uinfcrossz_array = cross.(eachrow(va_array), eachrow(z_airf_array)) + relative_velocity_array .= va_array .+ induced_velocity_all + for i in 1:n_panels + relative_velocity_crossz[i, :] .= cross3( + view(relative_velocity_array, i, :), + view(z_airf_array, i, :) + ) + Uinfcrossz_array[i, :] .= cross3( + view(va_array, i, :), + view(z_airf_array, i, :) + ) + end - v_normal_array = vec(sum(x_airf_array .* relative_velocity_array, dims=2)) - v_tangential_array = vec(sum(y_airf_array .* relative_velocity_array, dims=2)) + for i in 1:n_panels + v_normal_array[i] = dot(view(x_airf_array, i, :), view(relative_velocity_array, i, :)) + v_tangential_array[i] = dot(view(y_airf_array, i, :), view(relative_velocity_array, i, :)) + end alpha_array .= atan.(v_normal_array, v_tangential_array) + + for i in 1:n_panels + @views Umag_array[i] = norm(relative_velocity_crossz[i, :]) + @views Umagw_array[i] = norm(Uinfcrossz_array[i, :]) + end - Umag_array .= norm.(relative_velocity_crossz) - Umagw_array = norm.(Uinfcrossz_array) - - cl_array = [calculate_cl(panel, alpha) for (panel, alpha) in zip(panels, alpha_array)] + for (i, (panel, alpha)) in enumerate(zip(panels, alpha_array)) + cl_array[i] = calculate_cl(panel, alpha) + end gamma_new .= 0.5 .* Umag_array.^2 ./ Umagw_array .* cl_array .* chord_array # Apply damping if needed @@ -229,18 +256,19 @@ function gamma_loop( damp, is_damping_applied = smooth_circulation(gamma, 0.1, 0.5) @debug "damp: $damp" else - damp = zeros(length(gamma)) + damp .= 0.0 is_damping_applied = false end - # Update gamma with relaxation and damping - gamma_new = (1 - relaxation_factor) .* gamma + + gamma_new .= (1 - relaxation_factor) .* gamma .+ relaxation_factor .* gamma_new .+ damp # Check convergence - reference_error = maximum(abs.(gamma_new)) + abs_gamma_new .= abs.(gamma_new) + reference_error = maximum(abs_gamma_new) reference_error = max(reference_error, solver.tol_reference_error) - error = maximum(abs.(gamma_new - gamma)) + abs_gamma_new .= abs.(gamma_new .- gamma) + error = maximum(abs_gamma_new) normalized_error = error / reference_error @debug "Iteration: $i, normalized_error: $normalized_error, is_damping_applied: $is_damping_applied" From 7b1c098a9bf4a4ae63afdcadbf1b0bfa3ac0c0ec Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Fri, 21 Feb 2025 09:48:06 -0600 Subject: [PATCH 02/24] add dat file --- data/HL5_ram_air_kite_foil.dat | 46 ++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 data/HL5_ram_air_kite_foil.dat diff --git a/data/HL5_ram_air_kite_foil.dat b/data/HL5_ram_air_kite_foil.dat new file mode 100644 index 00000000..09eceaee --- /dev/null +++ b/data/HL5_ram_air_kite_foil.dat @@ -0,0 +1,46 @@ +A + 1.0000000 0.0000479 + 0.9666341 0.0017363 + 0.9038341 0.0066154 + 0.8342935 0.0136262 + 0.7624304 0.0214202 + 0.6907649 0.0299758 + 0.6203655 0.0397261 + 0.5500472 0.0504400 + 0.4788028 0.0614558 + 0.4074159 0.0732084 + 0.3368400 0.0848204 + 0.2692757 0.0941998 + 0.2083741 0.0988242 + 0.1573868 0.0983570 + 0.1169140 0.0924067 + 0.0828168 0.0826568 + 0.0560463 0.0704667 + 0.0361758 0.0577955 + 0.0215498 0.0452424 + 0.0115826 0.0329061 + 0.0055625 0.0219130 + 0.0019945 0.0124293 + 0.0002467 0.0035649 + 0.0004755 -0.0045913 + 0.0027675 -0.0130849 + 0.0077900 -0.0220373 + 0.0162007 -0.0306320 + 0.0285895 -0.0396469 + 0.0466671 -0.0494062 + 0.0728758 -0.0593274 + 0.1094090 -0.0681156 + 0.1556860 -0.0747603 + 0.2093904 -0.0783204 + 0.2712652 -0.0798598 + 0.3359510 -0.0788043 + 0.4038697 -0.0747532 + 0.4742506 -0.0687364 + 0.5453696 -0.0615892 + 0.6158254 -0.0528273 + 0.6875958 -0.0427482 + 0.7594140 -0.0320344 + 0.8304601 -0.0220930 + 0.9008044 -0.0130373 + 0.9663063 -0.0044401 + 1.0000000 0.0000479 From 605f8f34af05accfc40e9142cdaa0fc2e9216a45 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Fri, 21 Feb 2025 21:52:45 -0600 Subject: [PATCH 03/24] add bench file --- .gitignore | 1 + src/filament.jl | 6 +-- test/bench.jl | 115 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 118 insertions(+), 4 deletions(-) create mode 100644 test/bench.jl diff --git a/.gitignore b/.gitignore index 7fc38c41..51764540 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ Manifest.toml venv results/TUDELFT_V3_LEI_KITE/polars/tutorial_testing_stall_model_n_panels_54_distribution_split_provided.pdf docs/build/ +*.bin diff --git a/src/filament.jl b/src/filament.jl index 1d6229b5..9223b67f 100644 --- a/src/filament.jl +++ b/src/filament.jl @@ -57,8 +57,8 @@ function velocity_3D_bound_vortex!( elseif norm(r1Xr0) / norm(r0) == 0 vel .= zeros(3) else - @info "inside core radius" - @info "distance from control point to filament: $(norm(r1Xr0) / norm(r0))" + @debug "inside core radius" + @debug "distance from control point to filament: $(norm(r1Xr0) / norm(r0))" # Project onto core radius cross3!(r2Xr0, r2, r0) @@ -71,8 +71,6 @@ function velocity_3D_bound_vortex!( vel_ind_proj .= (gamma / (4π)) .* r1_projXr2_proj ./ (norm(r1_projXr2_proj)^2) .* dot(r0, r1_proj/norm(r1_proj) .- r2_proj/norm(r2_proj)) - @show vel_ind_proj r1_projXr2_proj - vel .= norm(r1Xr0) ./ (norm(r0) * epsilon) .* vel_ind_proj end nothing diff --git a/test/bench.jl b/test/bench.jl new file mode 100644 index 00000000..392bcc3b --- /dev/null +++ b/test/bench.jl @@ -0,0 +1,115 @@ +using BenchmarkTools +using VortexStepMethod +using VortexStepMethod: calculate_AIC_matrices, gamma_loop, calculate_results, + update_effective_angle_of_attack_if_VSM, calculate_projected_area, + calculate_cl, calculate_cd_cm, + calculate_velocity_induced_single_ring_semiinfinite, + calculate_velocity_induced_bound_2D, + velocity_3D_bound_vortex!, + velocity_3D_trailing_vortex!, + velocity_3D_trailing_vortex_semiinfinite!, + Panel +using Test +using LinearAlgebra + +@testset "Function Allocation Tests" begin + # 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) + + # Create test panels + panels = [] + 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 + "inviscid") + add_section!(wing, + [0.0, -span/2, 0.0], # Right tip LE + [chord, -span/2, 0.0], # Right tip TE + "inviscid") + + wing_aero = WingAerodynamics([wing]) + + vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a + set_va!(wing_aero, (vel_app, 0.0)) # Second parameter is yaw rate + + # Initialize solvers for both LLT and VSM methods + solver = Solver() + + # Pre-allocate arrays + gamma = rand(n_panels) + gamma_new = similar(gamma) + AIC_x = rand(n_panels, n_panels) + AIC_y = similar(AIC_x) + AIC_z = similar(AIC_x) + v_ind = zeros(3) + point = rand(3) + va_norm_array = ones(n_panels) + va_unit_array = ones(n_panels, 3) + + models = ["VSM", "LLT"] + core_radius_fractions = [0.001, 10.0] + + @testset "AIC Matrix Calculation" begin + for model in models + for frac in core_radius_fractions + result = @benchmark calculate_AIC_matrices($wing_aero, $model, $frac, $va_norm_array, $va_unit_array) + @test result.allocs ≤ 100 # Allow some allocations for matrix setup + end + end + end + + @testset "Gamma Loop" begin + result = @benchmark gamma_loop($solver, $wing, $gamma_new, $AIC_x, $AIC_y, $AIC_z) + @test result.allocs ≤ 50 # Main iteration should be mostly allocation-free + end + + @testset "Results Calculation" begin + result = @benchmark calculate_results($wing, $gamma) + @test result.allocs ≤ 20 # Allow minimal allocations for results + end + + @testset "Angle of Attack Update" begin + result = @benchmark update_effective_angle_of_attack_if_VSM($wing, $gamma) + @test result.allocs == 0 # Should be allocation-free + end + + @testset "Area Calculations" begin + result = @benchmark calculate_projected_area($wing) + @test result.allocs ≤ 10 # Geometric calculations may need some allocations + end + + @testset "Aerodynamic Coefficients" begin + panel = panels[1] + alpha = 0.1 + + @test (@ballocated calculate_cl($panel, $alpha)) == 0 + @test (@ballocated calculate_cd_cm($panel, $alpha)) == 0 + end + + @testset "Induced Velocity Calculations" begin + # Test single ring velocity calculation + @test (@ballocated calculate_velocity_induced_single_ring_semiinfinite( + $point, $panels[1], $gamma[1])) == 0 + + # Test 2D bound vortex + @test (@ballocated calculate_velocity_induced_bound_2D( + $point, $panels[1], $gamma[1])) == 0 + + # Test 3D velocity components + @test (@ballocated velocity_3D_bound_vortex!( + $v_ind, $point, $panels[1], $gamma[1])) == 0 + + @test (@ballocated velocity_3D_trailing_vortex!( + $v_ind, $point, $panels[1], $gamma[1])) == 0 + + @test (@ballocated velocity_3D_trailing_vortex_semiinfinite!( + $v_ind, $point, $panels[1], $gamma[1])) == 0 + end +end \ No newline at end of file From aaca68ffa41908db08001242686da99876c2da5e Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 22 Feb 2025 08:01:22 -0600 Subject: [PATCH 04/24] improve test --- test/bench.jl | 80 ++++++++++++++++++++++++++------------------------- 1 file changed, 41 insertions(+), 39 deletions(-) diff --git a/test/bench.jl b/test/bench.jl index 392bcc3b..43c6396d 100644 --- a/test/bench.jl +++ b/test/bench.jl @@ -59,57 +59,59 @@ using LinearAlgebra @testset "AIC Matrix Calculation" begin for model in models for frac in core_radius_fractions - result = @benchmark calculate_AIC_matrices($wing_aero, $model, $frac, $va_norm_array, $va_unit_array) - @test result.allocs ≤ 100 # Allow some allocations for matrix setup + @testset "Model $model Core Radius Fraction $frac" begin + result = @benchmark calculate_AIC_matrices($wing_aero, $model, $frac, $va_norm_array, $va_unit_array) + @test result.allocs ≤ 100 # Allow some allocations for matrix setup + end end end end - @testset "Gamma Loop" begin - result = @benchmark gamma_loop($solver, $wing, $gamma_new, $AIC_x, $AIC_y, $AIC_z) - @test result.allocs ≤ 50 # Main iteration should be mostly allocation-free - end + # @testset "Gamma Loop" begin + # result = @benchmark gamma_loop($solver, $wing, $gamma_new, $AIC_x, $AIC_y, $AIC_z) + # @test result.allocs ≤ 50 # Main iteration should be mostly allocation-free + # end - @testset "Results Calculation" begin - result = @benchmark calculate_results($wing, $gamma) - @test result.allocs ≤ 20 # Allow minimal allocations for results - end + # @testset "Results Calculation" begin + # result = @benchmark calculate_results($wing, $gamma) + # @test result.allocs ≤ 20 # Allow minimal allocations for results + # end - @testset "Angle of Attack Update" begin - result = @benchmark update_effective_angle_of_attack_if_VSM($wing, $gamma) - @test result.allocs == 0 # Should be allocation-free - end + # @testset "Angle of Attack Update" begin + # result = @benchmark update_effective_angle_of_attack_if_VSM($wing, $gamma) + # @test result.allocs == 0 # Should be allocation-free + # end - @testset "Area Calculations" begin - result = @benchmark calculate_projected_area($wing) - @test result.allocs ≤ 10 # Geometric calculations may need some allocations - end + # @testset "Area Calculations" begin + # result = @benchmark calculate_projected_area($wing) + # @test result.allocs ≤ 10 # Geometric calculations may need some allocations + # end - @testset "Aerodynamic Coefficients" begin - panel = panels[1] - alpha = 0.1 + # @testset "Aerodynamic Coefficients" begin + # panel = panels[1] + # alpha = 0.1 - @test (@ballocated calculate_cl($panel, $alpha)) == 0 - @test (@ballocated calculate_cd_cm($panel, $alpha)) == 0 - end + # @test (@ballocated calculate_cl($panel, $alpha)) == 0 + # @test (@ballocated calculate_cd_cm($panel, $alpha)) == 0 + # end - @testset "Induced Velocity Calculations" begin - # Test single ring velocity calculation - @test (@ballocated calculate_velocity_induced_single_ring_semiinfinite( - $point, $panels[1], $gamma[1])) == 0 + # @testset "Induced Velocity Calculations" begin + # # Test single ring velocity calculation + # @test (@ballocated calculate_velocity_induced_single_ring_semiinfinite( + # $point, $panels[1], $gamma[1])) == 0 - # Test 2D bound vortex - @test (@ballocated calculate_velocity_induced_bound_2D( - $point, $panels[1], $gamma[1])) == 0 + # # Test 2D bound vortex + # @test (@ballocated calculate_velocity_induced_bound_2D( + # $point, $panels[1], $gamma[1])) == 0 - # Test 3D velocity components - @test (@ballocated velocity_3D_bound_vortex!( - $v_ind, $point, $panels[1], $gamma[1])) == 0 + # # Test 3D velocity components + # @test (@ballocated velocity_3D_bound_vortex!( + # $v_ind, $point, $panels[1], $gamma[1])) == 0 - @test (@ballocated velocity_3D_trailing_vortex!( - $v_ind, $point, $panels[1], $gamma[1])) == 0 + # @test (@ballocated velocity_3D_trailing_vortex!( + # $v_ind, $point, $panels[1], $gamma[1])) == 0 - @test (@ballocated velocity_3D_trailing_vortex_semiinfinite!( - $v_ind, $point, $panels[1], $gamma[1])) == 0 - end + # @test (@ballocated velocity_3D_trailing_vortex_semiinfinite!( + # $v_ind, $point, $panels[1], $gamma[1])) == 0 + # end end \ No newline at end of file From 650bf72d519f64ddf2b2e3a0e6d8d435266c8180 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 22 Feb 2025 10:20:22 -0600 Subject: [PATCH 05/24] rename Umag to v_a --- data/HL5_ram_air_kite_foil.dat | 46 ---------------------------------- examples/rectangular_wing.jl | 8 +++--- src/solver.jl | 8 +++--- 3 files changed, 9 insertions(+), 53 deletions(-) delete mode 100644 data/HL5_ram_air_kite_foil.dat diff --git a/data/HL5_ram_air_kite_foil.dat b/data/HL5_ram_air_kite_foil.dat deleted file mode 100644 index 09eceaee..00000000 --- a/data/HL5_ram_air_kite_foil.dat +++ /dev/null @@ -1,46 +0,0 @@ -A - 1.0000000 0.0000479 - 0.9666341 0.0017363 - 0.9038341 0.0066154 - 0.8342935 0.0136262 - 0.7624304 0.0214202 - 0.6907649 0.0299758 - 0.6203655 0.0397261 - 0.5500472 0.0504400 - 0.4788028 0.0614558 - 0.4074159 0.0732084 - 0.3368400 0.0848204 - 0.2692757 0.0941998 - 0.2083741 0.0988242 - 0.1573868 0.0983570 - 0.1169140 0.0924067 - 0.0828168 0.0826568 - 0.0560463 0.0704667 - 0.0361758 0.0577955 - 0.0215498 0.0452424 - 0.0115826 0.0329061 - 0.0055625 0.0219130 - 0.0019945 0.0124293 - 0.0002467 0.0035649 - 0.0004755 -0.0045913 - 0.0027675 -0.0130849 - 0.0077900 -0.0220373 - 0.0162007 -0.0306320 - 0.0285895 -0.0396469 - 0.0466671 -0.0494062 - 0.0728758 -0.0593274 - 0.1094090 -0.0681156 - 0.1556860 -0.0747603 - 0.2093904 -0.0783204 - 0.2712652 -0.0798598 - 0.3359510 -0.0788043 - 0.4038697 -0.0747532 - 0.4742506 -0.0687364 - 0.5453696 -0.0615892 - 0.6158254 -0.0528273 - 0.6875958 -0.0427482 - 0.7594140 -0.0320344 - 0.8304601 -0.0220930 - 0.9008044 -0.0130373 - 0.9663063 -0.0044401 - 1.0000000 0.0000479 diff --git a/examples/rectangular_wing.jl b/examples/rectangular_wing.jl index 134cd9d4..477ce826 100644 --- a/examples/rectangular_wing.jl +++ b/examples/rectangular_wing.jl @@ -2,6 +2,8 @@ using LinearAlgebra using ControlPlots using VortexStepMethod +plot = false + # Step 1: Define wing parameters n_panels = 20 # Number of panels span = 20.0 # Wing span [m] @@ -53,7 +55,7 @@ println("CD = $(round(results_vsm["cd"], digits=4))") println("Projected area = $(round(results_vsm["projected_area"], digits=4)) m²") # Step 6: Plot geometry -plot_geometry( +plot && plot_geometry( wa, "Rectangular_wing_geometry"; data_type=".pdf", @@ -66,7 +68,7 @@ nothing # Step 7: Plot spanwise distributions y_coordinates = [panel.aerodynamic_center[2] for panel in wa.panels] -plot_distribution( +plot && plot_distribution( [y_coordinates, y_coordinates], [results_vsm, results_llt], ["VSM", "LLT"], @@ -75,7 +77,7 @@ plot_distribution( # Step 8: Plot polar curves angle_range = range(0, 20, 20) -plot_polars( +plot && plot_polars( [llt_solver, vsm_solver], [wa, wa], ["LLT", "VSM"]; diff --git a/src/solver.jl b/src/solver.jl index e98db9c1..30e43556 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -199,8 +199,8 @@ function gamma_loop( converged = false n_panels = wing_aero.n_panels alpha_array = wing_aero.alpha_array - Umag_array = wing_aero.Umag_array - Umagw_array = similar(Umag_array) + v_a_array = wing_aero.v_a_array + Umagw_array = similar(v_a_array) gamma = copy(gamma_new) abs_gamma_new = copy(gamma_new) @@ -242,14 +242,14 @@ function gamma_loop( alpha_array .= atan.(v_normal_array, v_tangential_array) for i in 1:n_panels - @views Umag_array[i] = norm(relative_velocity_crossz[i, :]) + @views v_a_array[i] = norm(relative_velocity_crossz[i, :]) @views Umagw_array[i] = norm(Uinfcrossz_array[i, :]) end for (i, (panel, alpha)) in enumerate(zip(panels, alpha_array)) cl_array[i] = calculate_cl(panel, alpha) end - gamma_new .= 0.5 .* Umag_array.^2 ./ Umagw_array .* cl_array .* chord_array + gamma_new .= 0.5 .* v_a_array.^2 ./ Umagw_array .* cl_array .* chord_array # Apply damping if needed if solver.is_with_artificial_damping From 0ab51528b2e02d38720d411db6fede8a641227e4 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 22 Feb 2025 11:37:07 -0600 Subject: [PATCH 06/24] cant remove calculate_velocity_induced_single_ring_semiinfinite! allocs --- src/filament.jl | 32 +++++++++++++------------- src/panel.jl | 32 +++++++++++++------------- src/solver.jl | 15 ++++--------- src/wing_aerodynamics.jl | 41 +++++++++++++++++++++------------- test/bench.jl | 4 ++-- test/test_wing_aerodynamics.jl | 5 +++-- 6 files changed, 66 insertions(+), 63 deletions(-) diff --git a/src/filament.jl b/src/filament.jl index 9223b67f..9d49d8ca 100644 --- a/src/filament.jl +++ b/src/filament.jl @@ -31,12 +31,12 @@ end Calculate induced velocity by a bound vortex filament at a point in space. """ function velocity_3D_bound_vortex!( - vel::VelVector, + vel, filament::BoundFilament, - XVP::PosVector, - gamma::Float64, - core_radius_fraction::Float64, - work_vectors::NTuple{10, Vector{Float64}} + XVP, + gamma, + core_radius_fraction, + work_vectors ) r1, r2, r1Xr2, r1Xr0, r2Xr0, r1r2norm, r1_proj, r2_proj, r1_projXr2_proj, vel_ind_proj = work_vectors r0 = filament.r0 @@ -93,12 +93,12 @@ Reference: Rick Damiani et al. "A vortex step method for nonlinear airfoil polar as implemented in KiteAeroDyn". """ function velocity_3D_trailing_vortex!( - vel::VelVector, + vel, filament::BoundFilament, - XVP::PosVector, - gamma::Float64, - Uinf::Float64, - work_vectors::NTuple{10,Vector{Float64}} + XVP, + gamma, + Uinf, + work_vectors ) r0, r1, r2, r_perp, r1Xr2, r1Xr0, r2Xr0, normr1r2 = work_vectors[1:8] r0 .= filament.x2 .- filament.x1 # Vortex filament @@ -160,13 +160,13 @@ end Calculate induced velocity by a semi-infinite trailing vortex filament. """ function velocity_3D_trailing_vortex_semiinfinite!( - vel::VelVector, + vel, filament::SemiInfiniteFilament, - Vf::VelVector, - XVP::PosVector, - GAMMA::Float64, - Uinf::Float64, - work_vectors::NTuple{10,Vector{Float64}} + Vf, + XVP, + GAMMA, + Uinf, + work_vectors ) r1, r_perp, r1XVf = work_vectors[1:3] GAMMA = -GAMMA * filament.filament_direction diff --git a/src/panel.jl b/src/panel.jl index f4694ad3..34946814 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -338,7 +338,9 @@ function calculate_filaments_for_plotting(panel::Panel) end """ - calculate_velocity_induced_single_ring_semiinfinite( + calculate_velocity_induced_single_ring_semiinfinite!( + velind::Matrix{Float64}, + tempvel::Vector{Float64}, panel::Panel, evaluation_point::Vector{Float64}, evaluation_point_on_bound::Bool, @@ -359,26 +361,25 @@ Calculate the velocity induced by a vortex ring at a control point. - `core_radius_fraction`: Vortex core radius as fraction of panel width # Returns -- `Vector{Float64}`: Induced velocity vector +- nothing """ -function calculate_velocity_induced_single_ring_semiinfinite( +function calculate_velocity_induced_single_ring_semiinfinite!( + velind, + tempvel, panel::Panel, - evaluation_point::PosVector, - evaluation_point_on_bound::Bool, - va_norm::Float64, - va_unit::VelVector, - gamma::Float64, - core_radius_fraction::Float64, - work_vectors::NTuple{10,Vector{Float64}} + evaluation_point, + evaluation_point_on_bound, + va_norm, + va_unit, + gamma, + core_radius_fraction, + work_vectors ) - velind = zeros(Float64, 3) - tempvel = zeros(3) - # Process each filament for (i, filament) in enumerate(panel.filaments) if i == 1 # bound filament if evaluation_point_on_bound - tempvel .= zeros(Float64, 3) + tempvel .= 0.0 else velocity_3D_bound_vortex!( tempvel, @@ -413,8 +414,7 @@ function calculate_velocity_induced_single_ring_semiinfinite( end velind .+= tempvel end - - return velind + return nothing end """ diff --git a/src/solver.jl b/src/solver.jl index 30e43556..f76a8399 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -89,7 +89,7 @@ function solve(solver::Solver, wing_aero::WingAerodynamics, gamma_distribution=n va_unit_array = va_array ./ va_norm_array # Calculate AIC matrices - AIC_x, AIC_y, AIC_z = calculate_AIC_matrices( + calculate_AIC_matrices!( wing_aero, solver.aerodynamic_model_type, solver.core_radius_fraction, @@ -116,9 +116,6 @@ function solve(solver::Solver, wing_aero::WingAerodynamics, gamma_distribution=n solver, wing_aero, gamma_initial, - AIC_x, - AIC_y, - AIC_z, va_array, chord_array, x_airf_array, @@ -134,9 +131,6 @@ function solve(solver::Solver, wing_aero::WingAerodynamics, gamma_distribution=n solver, wing_aero, gamma_initial, - AIC_x, - AIC_y, - AIC_z, va_array, chord_array, x_airf_array, @@ -185,9 +179,6 @@ function gamma_loop( solver::Solver, wing_aero::WingAerodynamics, gamma_new::Vector{Float64}, - AIC_x::Matrix{Float64}, - AIC_y::Matrix{Float64}, - AIC_z::Matrix{Float64}, va_array::Matrix{Float64}, chord_array::Vector{Float64}, x_airf_array::Matrix{Float64}, @@ -204,7 +195,7 @@ function gamma_loop( gamma = copy(gamma_new) abs_gamma_new = copy(gamma_new) - induced_velocity_all = zeros(size(AIC_x, 1), 3) + induced_velocity_all = zeros(n_panels, 3) relative_velocity_array = similar(va_array) relative_velocity_crossz = similar(relative_velocity_array) Uinfcrossz_array = similar(va_array) @@ -212,6 +203,8 @@ function gamma_loop( damp = zeros(length(gamma)) v_normal_array = zeros(n_panels) v_tangential_array = zeros(n_panels) + + AIC_x, AIC_y, AIC_z = @views wing_aero.AIC[1, :, :], wing_aero.AIC[2, :, :], wing_aero.AIC[3, :, :] iters = 0 for i in 1:solver.max_iterations diff --git a/src/wing_aerodynamics.jl b/src/wing_aerodynamics.jl index d646186c..a0447ca5 100644 --- a/src/wing_aerodynamics.jl +++ b/src/wing_aerodynamics.jl @@ -16,7 +16,8 @@ mutable struct WingAerodynamics alpha_array::Vector{Float64} v_a_array::Vector{Float64} - work_vectors::NTuple{10,Vector{Float64}} + work_vectors::NTuple{10,MVec3} + AIC::Array{Float64, 3} function WingAerodynamics( wings::Vector{T}; @@ -60,7 +61,8 @@ mutable struct WingAerodynamics alpha_array = zeros(n_panels) v_a_array = zeros(n_panels) work_vectors = ntuple(_ -> Vector{Float64}(undef, 3), 10) - + AIC = zeros(3, n_panels, n_panels) + new( panels, n_panels, @@ -72,7 +74,8 @@ mutable struct WingAerodynamics stall_angle_list, alpha_array, v_a_array, - work_vectors + work_vectors, + AIC ) end end @@ -206,7 +209,7 @@ function calculate_panel_properties(section_list::Vector{Section}, n_panels::Int end """ - calculate_AIC_matrices(wa::WingAerodynamics, model::String, + calculate_AIC_matrices!(wa::WingAerodynamics, model::String, core_radius_fraction::Float64, va_norm_array::Vector{Float64}, va_unit_array::Matrix{Float64}) @@ -216,35 +219,41 @@ Calculate Aerodynamic Influence Coefficient matrices. Returns: Tuple of (AIC_x, AIC_y, AIC_z) matrices """ -function calculate_AIC_matrices(wa::WingAerodynamics, model::String, - core_radius_fraction::Float64, - va_norm_array::Vector{Float64}, - va_unit_array::Matrix{Float64}) +function calculate_AIC_matrices!(wa::WingAerodynamics, model, + core_radius_fraction, + va_norm_array, + va_unit_array) model in ["VSM", "LLT"] || throw(ArgumentError("Model must be VSM or LLT")) - # Determine evaluation point based on model evaluation_point = model == "VSM" ? :control_point : :aerodynamic_center evaluation_point_on_bound = model == "LLT" # Initialize AIC matrices - AIC = zeros(3, wa.n_panels, wa.n_panels) + AIC = @views wa.AIC + velocity_induced, tempvel, va_unit = zeros(MVec3), zeros(MVec3), zeros(MVec3) # Calculate influence coefficients for icp in 1:wa.n_panels ep = getproperty(wa.panels[icp], evaluation_point) for jring in 1:wa.n_panels - velocity_induced = calculate_velocity_induced_single_ring_semiinfinite( + velocity_induced .= 0.0 + tempvel .= 0.0 + va_unit .= @views va_unit_array[jring, :] + wa.panels[jring] + calculate_velocity_induced_single_ring_semiinfinite!( + velocity_induced, + tempvel, wa.panels[jring], ep, evaluation_point_on_bound, va_norm_array[jring], - va_unit_array[jring, :], + va_unit, 1.0, core_radius_fraction, wa.work_vectors ) - AIC[:, icp, jring] = velocity_induced + AIC[:, icp, jring] .= velocity_induced # Subtract 2D induced velocity for VSM if icp == jring && model == "VSM" @@ -253,8 +262,7 @@ function calculate_AIC_matrices(wa::WingAerodynamics, model::String, end end end - - return AIC[1, :, :], AIC[2, :, :], AIC[3, :, :] + return nothing end """ @@ -352,9 +360,10 @@ function update_effective_angle_of_attack_if_VSM(wa::WingAerodynamics, va_unit_array::Matrix{Float64}) # Calculate AIC matrices at aerodynamic center using LLT method - AIC_x, AIC_y, AIC_z = calculate_AIC_matrices( + calculate_AIC_matrices!( wa, "LLT", core_radius_fraction, va_norm_array, va_unit_array ) + AIC_x, AIC_y, AIC_z = @views wa.AIC[1, :, :], wa.AIC[2, :, :], wa.AIC[3, :, :] # Calculate induced velocities induced_velocity = [ diff --git a/test/bench.jl b/test/bench.jl index 43c6396d..a9a871a6 100644 --- a/test/bench.jl +++ b/test/bench.jl @@ -1,6 +1,6 @@ using BenchmarkTools using VortexStepMethod -using VortexStepMethod: calculate_AIC_matrices, gamma_loop, calculate_results, +using VortexStepMethod: calculate_AIC_matrices!, gamma_loop, calculate_results, update_effective_angle_of_attack_if_VSM, calculate_projected_area, calculate_cl, calculate_cd_cm, calculate_velocity_induced_single_ring_semiinfinite, @@ -60,7 +60,7 @@ using LinearAlgebra for model in models for frac in core_radius_fractions @testset "Model $model Core Radius Fraction $frac" begin - result = @benchmark calculate_AIC_matrices($wing_aero, $model, $frac, $va_norm_array, $va_unit_array) + result = @benchmark calculate_AIC_matrices!($wing_aero, $model, $frac, $va_norm_array, $va_unit_array) @test result.allocs ≤ 100 # Allow some allocations for matrix setup end end diff --git a/test/test_wing_aerodynamics.jl b/test/test_wing_aerodynamics.jl index 8014ca6b..6a36497a 100644 --- a/test/test_wing_aerodynamics.jl +++ b/test/test_wing_aerodynamics.jl @@ -1,6 +1,6 @@ # using VortexStepMethod: Wing, WingAerodynamics, BoundFilament, SemiInfiniteFilament, add_section!, set_va!, solve, calculate_cl using VortexStepMethod -using VortexStepMethod: calculate_cl, calculate_cd_cm, calculate_projected_area, calculate_AIC_matrices +using VortexStepMethod: calculate_cl, calculate_cd_cm, calculate_projected_area, calculate_AIC_matrices! using LinearAlgebra using Test using Logging @@ -159,13 +159,14 @@ end # Calculate new matrices va_norm_array = fill(norm(Uinf), length(coord)) va_unit_array = repeat(reshape(Uinf ./ norm(Uinf), 1, 3), length(coord)) - AIC_x, AIC_y, AIC_z = calculate_AIC_matrices( + calculate_AIC_matrices!( wing_aero, "LLT", core_radius_fraction, va_norm_array, va_unit_array ) + AIC_x, AIC_y, AIC_z = @views wing_aero.AIC[1, :, :], wing_aero.AIC[2, :, :], wing_aero.AIC[3, :, :] # Compare matrices @test isapprox(MatrixU, AIC_x, atol=1e-5) From 7e2bdbcc1d0b0de19e349d90a18b846b2a400dfa Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 22 Feb 2025 17:48:15 -0600 Subject: [PATCH 07/24] two non-allocating functions in code_bench --- src/VortexStepMethod.jl | 1 + src/filament.jl | 2 +- src/panel.jl | 54 ++++++++++++++++++++++------------------ src/wing_aerodynamics.jl | 13 +++++----- 4 files changed, 39 insertions(+), 31 deletions(-) diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 33d47936..e9076c6d 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -14,6 +14,7 @@ using Interpolations: linear_interpolation, Line, Extrapolation using Serialization using Distributed using SharedArrays +using BenchmarkTools # Export public interface export Wing, Section, KiteWing diff --git a/src/filament.jl b/src/filament.jl index 9d49d8ca..3e9f0afe 100644 --- a/src/filament.jl +++ b/src/filament.jl @@ -92,7 +92,7 @@ Calculate induced velocity by a trailing vortex filament. Reference: Rick Damiani et al. "A vortex step method for nonlinear airfoil polar data as implemented in KiteAeroDyn". """ -function velocity_3D_trailing_vortex!( +@inline function velocity_3D_trailing_vortex!( vel, filament::BoundFilament, XVP, diff --git a/src/panel.jl b/src/panel.jl index 34946814..f6e2e4d1 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -363,10 +363,10 @@ Calculate the velocity induced by a vortex ring at a control point. # Returns - nothing """ -function calculate_velocity_induced_single_ring_semiinfinite!( +@inline function calculate_velocity_induced_single_ring_semiinfinite!( velind, tempvel, - panel::Panel, + filaments, evaluation_point, evaluation_point_on_bound, va_norm, @@ -375,42 +375,44 @@ function calculate_velocity_induced_single_ring_semiinfinite!( core_radius_fraction, work_vectors ) + velind .= 0.0 + tempvel .= 0.0 # Process each filament - for (i, filament) in enumerate(panel.filaments) + @inbounds for i in eachindex(filaments) if i == 1 # bound filament if evaluation_point_on_bound tempvel .= 0.0 else velocity_3D_bound_vortex!( tempvel, - filament, + filaments[i], evaluation_point, gamma, core_radius_fraction, work_vectors ) end - elseif i ∈ (2, 3) # trailing filaments + elseif i == 2 || i == 3 # trailing filaments velocity_3D_trailing_vortex!( tempvel, - filament, - evaluation_point, - gamma, - va_norm, - work_vectors - ) - elseif i ∈ (4, 5) # semi-infinite trailing filaments - velocity_3D_trailing_vortex_semiinfinite!( - tempvel, - filament, - va_unit, + filaments[i], evaluation_point, gamma, va_norm, work_vectors ) + elseif i == 4 || i == 5 # semi-infinite trailing filaments + # velocity_3D_trailing_vortex_semiinfinite!( + # tempvel, + # filaments[i], + # va_unit, + # evaluation_point, + # gamma, + # va_norm, + # work_vectors + # ) else - tempvel .= zeros(Float64, 3) + tempvel .= 0.0 end velind .+= tempvel end @@ -418,7 +420,7 @@ function calculate_velocity_induced_single_ring_semiinfinite!( end """ - calculate_velocity_induced_bound_2D(panel::Panel, evaluation_point::Vector{Float64}) + calculate_velocity_induced_bound_2D!(panel::Panel, evaluation_point::Vector{Float64}) Calculate velocity induced by bound vortex filaments at the control point. Only needed for VSM, as LLT bound and filament align, thus no induced velocity. @@ -430,19 +432,23 @@ Only needed for VSM, as LLT bound and filament align, thus no induced velocity. # Returns - `Vector{Float64}`: Induced velocity at the control point """ -function calculate_velocity_induced_bound_2D( +function calculate_velocity_induced_bound_2D!( + U_2D, panel::Panel, - evaluation_point::PosVector + evaluation_point, + work_vectors ) + r3, r0, cross_, cross_square = work_vectors # r3 perpendicular to the bound vortex - r3 = evaluation_point - (panel.bound_point_1 + panel.bound_point_2) / 2 + r3 .= evaluation_point .- (panel.bound_point_1 .+ panel.bound_point_2) ./ 2 # r0 is the direction of the bound vortex - r0 = panel.bound_point_1 - panel.bound_point_2 + r0 .= panel.bound_point_1 .- panel.bound_point_2 # Calculate cross product - cross_ = cross(r0, r3) + cross3!(cross_, r0, r3) # Calculate induced velocity - return (cross_ / sum(cross_.^2) / 2π) * norm(r0) + U_2D .= (cross_ ./ sum(cross_square) ./ 2π) .* norm(r0) + return nothing end \ No newline at end of file diff --git a/src/wing_aerodynamics.jl b/src/wing_aerodynamics.jl index a0447ca5..c95ac27f 100644 --- a/src/wing_aerodynamics.jl +++ b/src/wing_aerodynamics.jl @@ -60,7 +60,7 @@ mutable struct WingAerodynamics alpha_array = zeros(n_panels) v_a_array = zeros(n_panels) - work_vectors = ntuple(_ -> Vector{Float64}(undef, 3), 10) + work_vectors = ntuple(_ -> zeros(MVec3), 10) AIC = zeros(3, n_panels, n_panels) new( @@ -230,7 +230,7 @@ function calculate_AIC_matrices!(wa::WingAerodynamics, model, # Initialize AIC matrices AIC = @views wa.AIC - velocity_induced, tempvel, va_unit = zeros(MVec3), zeros(MVec3), zeros(MVec3) + velocity_induced, tempvel, va_unit, U_2D = zeros(MVec3), zeros(MVec3), zeros(MVec3), zeros(MVec3) # Calculate influence coefficients for icp in 1:wa.n_panels @@ -239,14 +239,15 @@ function calculate_AIC_matrices!(wa::WingAerodynamics, model, velocity_induced .= 0.0 tempvel .= 0.0 va_unit .= @views va_unit_array[jring, :] - wa.panels[jring] + filaments = wa.panels[jring].filaments + va_norm = va_norm_array[jring] calculate_velocity_induced_single_ring_semiinfinite!( velocity_induced, tempvel, - wa.panels[jring], + filaments, ep, evaluation_point_on_bound, - va_norm_array[jring], + va_norm, va_unit, 1.0, core_radius_fraction, @@ -257,7 +258,7 @@ function calculate_AIC_matrices!(wa::WingAerodynamics, model, # Subtract 2D induced velocity for VSM if icp == jring && model == "VSM" - U_2D = calculate_velocity_induced_bound_2D(wa.panels[jring], ep) + calculate_velocity_induced_bound_2D!(U_2D, wa.panels[jring], ep, work_vectors) AIC[:, icp, jring] .-= U_2D end end From 74f68d4b508843122139302ac42a19d7d944a690 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 22 Feb 2025 17:49:02 -0600 Subject: [PATCH 08/24] add simple script for alloc testing --- examples/code_bench.jl | 75 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100644 examples/code_bench.jl diff --git a/examples/code_bench.jl b/examples/code_bench.jl new file mode 100644 index 00000000..d257e934 --- /dev/null +++ b/examples/code_bench.jl @@ -0,0 +1,75 @@ +using Revise +using LinearAlgebra +using ControlPlots +using VortexStepMethod + +# 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 = WingAerodynamics([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 +vsm_solver = Solver(aerodynamic_model_type="VSM") + +# Benchmark setup +velocity_induced = zeros(3) +tempvel = zeros(3) +panel = wa.panels[1] +ep = [0.25, 9.5, 0.0] +evaluation_point_on_bound = true +va_norm = 20.0 +va_unit = [0.8660254037844387, 0.0, 0.4999999999999999] +core_radius_fraction = 1.0e-20 +gamma = 1.0 + +# Create work vectors tuple +work_vectors = ntuple(_ -> zeros(3), 10) + +using BenchmarkTools +@btime VortexStepMethod.calculate_velocity_induced_single_ring_semiinfinite!( + $velocity_induced, + $tempvel, + $panel.filaments, + $ep, + $evaluation_point_on_bound, + $va_norm, + $va_unit, + $gamma, + $core_radius_fraction, + $work_vectors +) + +U_2D = zeros(3) + +@btime VortexStepMethod.calculate_velocity_induced_bound_2D!( + $U_2D, + $panel, + $ep, + $work_vectors +) + +nothing From 9bdea5ef16b5befca2af92a1f47d50da952ccb64 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 22 Feb 2025 17:58:07 -0600 Subject: [PATCH 09/24] fix calculations --- src/panel.jl | 20 ++++++++++---------- src/wing_aerodynamics.jl | 2 +- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/panel.jl b/src/panel.jl index f6e2e4d1..d2c92755 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -402,15 +402,15 @@ Calculate the velocity induced by a vortex ring at a control point. work_vectors ) elseif i == 4 || i == 5 # semi-infinite trailing filaments - # velocity_3D_trailing_vortex_semiinfinite!( - # tempvel, - # filaments[i], - # va_unit, - # evaluation_point, - # gamma, - # va_norm, - # work_vectors - # ) + velocity_3D_trailing_vortex_semiinfinite!( + tempvel, + filaments[i], + va_unit, + evaluation_point, + gamma, + va_norm, + work_vectors + ) else tempvel .= 0.0 end @@ -449,6 +449,6 @@ function calculate_velocity_induced_bound_2D!( cross3!(cross_, r0, r3) # Calculate induced velocity - U_2D .= (cross_ ./ sum(cross_square) ./ 2π) .* norm(r0) + U_2D .= -(cross_ ./ sum(cross_square) ./ 2π) .* norm(r0) return nothing end \ No newline at end of file diff --git a/src/wing_aerodynamics.jl b/src/wing_aerodynamics.jl index c95ac27f..2e33990b 100644 --- a/src/wing_aerodynamics.jl +++ b/src/wing_aerodynamics.jl @@ -258,7 +258,7 @@ function calculate_AIC_matrices!(wa::WingAerodynamics, model, # Subtract 2D induced velocity for VSM if icp == jring && model == "VSM" - calculate_velocity_induced_bound_2D!(U_2D, wa.panels[jring], ep, work_vectors) + calculate_velocity_induced_bound_2D!(U_2D, wa.panels[jring], ep, wa.work_vectors) AIC[:, icp, jring] .-= U_2D end end From 7a0a093f24276d38cbd4dbd16c88cb2aa4239d34 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 22 Feb 2025 19:36:53 -0600 Subject: [PATCH 10/24] three non-allocating functions --- examples/code_bench.jl | 48 +++++++++++++++++++++++++++++----------- src/panel.jl | 20 ++++++++--------- src/wing_aerodynamics.jl | 8 ++----- 3 files changed, 47 insertions(+), 29 deletions(-) diff --git a/examples/code_bench.jl b/examples/code_bench.jl index d257e934..359c4845 100644 --- a/examples/code_bench.jl +++ b/examples/code_bench.jl @@ -2,6 +2,8 @@ using Revise using LinearAlgebra using ControlPlots using VortexStepMethod +using StaticArrays +using BenchmarkTools # Step 1: Define wing parameters n_panels = 20 # Number of panels @@ -36,20 +38,19 @@ set_va!(wa, (vel_app, 0.0)) # Second parameter is yaw rate vsm_solver = Solver(aerodynamic_model_type="VSM") # Benchmark setup -velocity_induced = zeros(3) -tempvel = zeros(3) +velocity_induced = @MVector zeros(3) # StaticArraysCore.MVector{3, Float64} +tempvel = @MVector zeros(3) # StaticArraysCore.MVector{3, Float64} panel = wa.panels[1] -ep = [0.25, 9.5, 0.0] -evaluation_point_on_bound = true -va_norm = 20.0 -va_unit = [0.8660254037844387, 0.0, 0.4999999999999999] -core_radius_fraction = 1.0e-20 -gamma = 1.0 +ep = @MVector [0.25, 9.5, 0.0] # StaticArraysCore.MVector{3, Float64} +evaluation_point_on_bound = true # Bool +va_norm = 20.0 # Float64 +va_unit = @MVector [0.8660254037844387, 0.0, 0.4999999999999999] # StaticArraysCore.MVector{3, Float64} +core_radius_fraction = 1.0e-20 # Float64 +gamma = 1.0 # Float64 -# Create work vectors tuple -work_vectors = ntuple(_ -> zeros(3), 10) +# Create work vectors tuple of MVector{3, Float64} +work_vectors = ntuple(_ -> @MVector(zeros(3)), 10) # NTuple{10, StaticArraysCore.MVector{3, Float64}} -using BenchmarkTools @btime VortexStepMethod.calculate_velocity_induced_single_ring_semiinfinite!( $velocity_induced, $tempvel, @@ -63,7 +64,20 @@ using BenchmarkTools $work_vectors ) -U_2D = zeros(3) +@btime VortexStepMethod.calculate_velocity_induced_single_ring_semiinfinite!( + velocity_induced, + tempvel, + panel.filaments, + ep, + evaluation_point_on_bound, + va_norm, + va_unit, + gamma, + core_radius_fraction, + work_vectors +) + +U_2D = @MVector zeros(3) # StaticArraysCore.MVector{3, Float64} @btime VortexStepMethod.calculate_velocity_induced_bound_2D!( $U_2D, @@ -72,4 +86,12 @@ U_2D = zeros(3) $work_vectors ) -nothing +# model = "VSM" +# va_norm_array = zeros(wa.n_panels) +# va_unit_array = zeros(wa.n_panels, 3) +# @time VortexStepMethod.calculate_AIC_matrices!(wa, model, +# core_radius_fraction, +# va_norm_array, +# va_unit_array) + +nothing \ No newline at end of file diff --git a/src/panel.jl b/src/panel.jl index d2c92755..40fe20fc 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -364,16 +364,16 @@ Calculate the velocity induced by a vortex ring at a control point. - nothing """ @inline function calculate_velocity_induced_single_ring_semiinfinite!( - velind, - tempvel, - filaments, - evaluation_point, - evaluation_point_on_bound, - va_norm, - va_unit, - gamma, - core_radius_fraction, - work_vectors + velind::MVector{3, Float64}, + tempvel::MVector{3, Float64}, + filaments::Vector{Union{VortexStepMethod.BoundFilament, VortexStepMethod.SemiInfiniteFilament}}, + evaluation_point::MVector{3, Float64}, + evaluation_point_on_bound::Bool, + va_norm::Float64, + va_unit::MVector{3, Float64}, + gamma::Float64, + core_radius_fraction::Float64, + work_vectors::NTuple{10, MVector{3, Float64}} ) velind .= 0.0 tempvel .= 0.0 diff --git a/src/wing_aerodynamics.jl b/src/wing_aerodynamics.jl index 2e33990b..d3a11d38 100644 --- a/src/wing_aerodynamics.jl +++ b/src/wing_aerodynamics.jl @@ -229,15 +229,12 @@ function calculate_AIC_matrices!(wa::WingAerodynamics, model, evaluation_point_on_bound = model == "LLT" # Initialize AIC matrices - AIC = @views wa.AIC velocity_induced, tempvel, va_unit, U_2D = zeros(MVec3), zeros(MVec3), zeros(MVec3), zeros(MVec3) # Calculate influence coefficients for icp in 1:wa.n_panels ep = getproperty(wa.panels[icp], evaluation_point) for jring in 1:wa.n_panels - velocity_induced .= 0.0 - tempvel .= 0.0 va_unit .= @views va_unit_array[jring, :] filaments = wa.panels[jring].filaments va_norm = va_norm_array[jring] @@ -253,13 +250,12 @@ function calculate_AIC_matrices!(wa::WingAerodynamics, model, core_radius_fraction, wa.work_vectors ) - - AIC[:, icp, jring] .= velocity_induced + wa.AIC[:, icp, jring] .= velocity_induced # Subtract 2D induced velocity for VSM if icp == jring && model == "VSM" calculate_velocity_induced_bound_2D!(U_2D, wa.panels[jring], ep, wa.work_vectors) - AIC[:, icp, jring] .-= U_2D + wa.AIC[:, icp, jring] .-= U_2D end end end From 786b2dc6b46492be685b3a866ba883b8f515926f Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 22 Feb 2025 21:31:02 -0600 Subject: [PATCH 11/24] make tests pass again --- src/filament.jl | 6 +++--- src/panel.jl | 13 +++++++------ test/test_wing_aerodynamics.jl | 3 ++- 3 files changed, 12 insertions(+), 10 deletions(-) diff --git a/src/filament.jl b/src/filament.jl index 3e9f0afe..a6297975 100644 --- a/src/filament.jl +++ b/src/filament.jl @@ -55,7 +55,7 @@ function velocity_3D_bound_vortex!( vel .= (gamma / (4π)) .* r1Xr2 ./ (norm(r1Xr2)^2) .* dot(r0, r1r2norm) elseif norm(r1Xr0) / norm(r0) == 0 - vel .= zeros(3) + vel .= 0.0 else @debug "inside core radius" @debug "distance from control point to filament: $(norm(r1Xr0) / norm(r0))" @@ -122,7 +122,7 @@ as implemented in KiteAeroDyn". vel .= (gamma / (4π)) .* r1Xr2 ./ (norm(r1Xr2)^2) .* dot(r0, normr1r2) elseif norm(r1Xr0) / norm(r0) == 0 - vel .= zeros(3) + vel .= 0.0 else # Project onto core radius r1_proj = dot(r1, r0) * r0 / (norm(r0)^2) + @@ -182,7 +182,7 @@ function velocity_3D_trailing_vortex_semiinfinite!( K = GAMMA / (4π) / norm(r1XVf)^2 * (1 + dot(r1, Vf) / norm(r1)) vel .= K .* r1XVf elseif norm(r1XVf) / norm(Vf) == 0 - vel .= zeros(3) + vel .= 0.0 else r1_proj = dot(r1, Vf) * Vf + epsilon * (r1/norm(r1) - Vf) / norm(r1/norm(r1) - Vf) diff --git a/src/panel.jl b/src/panel.jl index 40fe20fc..80a1ce8e 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -364,16 +364,16 @@ Calculate the velocity induced by a vortex ring at a control point. - nothing """ @inline function calculate_velocity_induced_single_ring_semiinfinite!( - velind::MVector{3, Float64}, - tempvel::MVector{3, Float64}, + velind::MVec3, + tempvel::MVec3, filaments::Vector{Union{VortexStepMethod.BoundFilament, VortexStepMethod.SemiInfiniteFilament}}, - evaluation_point::MVector{3, Float64}, + evaluation_point::MVec3, evaluation_point_on_bound::Bool, va_norm::Float64, - va_unit::MVector{3, Float64}, + va_unit::MVec3, gamma::Float64, core_radius_fraction::Float64, - work_vectors::NTuple{10, MVector{3, Float64}} + work_vectors::NTuple{10, MVec3} ) velind .= 0.0 tempvel .= 0.0 @@ -449,6 +449,7 @@ function calculate_velocity_induced_bound_2D!( cross3!(cross_, r0, r3) # Calculate induced velocity - U_2D .= -(cross_ ./ sum(cross_square) ./ 2π) .* norm(r0) + cross_square .= cross_.^2 + U_2D .= (cross_ ./ sum(cross_square) ./ 2π) .* norm(r0) return nothing end \ No newline at end of file diff --git a/test/test_wing_aerodynamics.jl b/test/test_wing_aerodynamics.jl index 6a36497a..c45181af 100644 --- a/test/test_wing_aerodynamics.jl +++ b/test/test_wing_aerodynamics.jl @@ -194,13 +194,14 @@ end # Calculate new matrices va_norm_array = fill(norm(Uinf), length(coord)) va_unit_array = repeat(reshape(Uinf ./ norm(Uinf), 1, 3), length(coord)) - AIC_x, AIC_y, AIC_z = calculate_AIC_matrices( + calculate_AIC_matrices!( wing_aero, "VSM", core_radius_fraction, va_norm_array, va_unit_array ) + AIC_x, AIC_y, AIC_z = wing_aero.AIC[1, :, :], wing_aero.AIC[2, :, :], wing_aero.AIC[3, :, :] # Compare matrices with higher precision for VSM @test isapprox(MatrixU, AIC_x, atol=1e-8) From 7317bc499baf9e10f0286fc1ef9156dc656cfc6e Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sun, 23 Feb 2025 09:23:28 -0600 Subject: [PATCH 12/24] non-allocating gamma loop --- examples/code_bench.jl | 15 +++++++------ src/panel.jl | 50 +++++++++++++++++++++++++----------------- src/solver.jl | 16 +++++++++----- src/wing_geometry.jl | 18 ++++++++++++--- 4 files changed, 63 insertions(+), 36 deletions(-) diff --git a/examples/code_bench.jl b/examples/code_bench.jl index 359c4845..cfe0f033 100644 --- a/examples/code_bench.jl +++ b/examples/code_bench.jl @@ -86,12 +86,13 @@ U_2D = @MVector zeros(3) # StaticArraysCore.MVector{3, Float64} $work_vectors ) -# model = "VSM" -# va_norm_array = zeros(wa.n_panels) -# va_unit_array = zeros(wa.n_panels, 3) -# @time VortexStepMethod.calculate_AIC_matrices!(wa, model, -# core_radius_fraction, -# va_norm_array, -# va_unit_array) +model = "VSM" +va_norm_array = zeros(wa.n_panels) +va_unit_array = zeros(wa.n_panels, 3) +@btime VortexStepMethod.calculate_AIC_matrices!( + $wa, model, + $core_radius_fraction, + $va_norm_array, + $va_unit_array) nothing \ No newline at end of file diff --git a/src/panel.jl b/src/panel.jl index 80a1ce8e..f23739de 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -24,7 +24,7 @@ Represents a panel in a vortex step method simulation. - `width::Float64`: Panel width - `filaments::Vector{BoundFilament}`: Panel filaments """ -mutable struct Panel{P} +mutable struct Panel TE_point_1::MVec3 LE_point_1::MVec3 TE_point_2::MVec3 @@ -36,7 +36,7 @@ mutable struct Panel{P} cl_coefficients::Union{Nothing,Vector{Float64}} cd_coefficients::Union{Nothing,Vector{Float64}} cm_coefficients::Union{Nothing,Vector{Float64}} - polar_data::P + polar_data::Union{Nothing, Tuple{Extrapolation}} aerodynamic_center::MVec3 control_point::MVec3 bound_point_1::MVec3 @@ -92,6 +92,22 @@ mutable struct Panel{P} 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 + ) + polar_data = (cl_interp, cd_interp, 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] @@ -113,7 +129,7 @@ mutable struct Panel{P} BoundFilament(TE_point_2, bound_point_2) ] - new{typeof(polar_data)}( + new( TE_point_1, LE_point_1, TE_point_2, LE_point_2, chord, nothing, corner_points, aero_model, cl_coeffs, cd_coeffs, cm_coeffs, polar_data, @@ -178,7 +194,7 @@ function compute_lei_coefficients(section_1::Section, section_2::Section) ) # Compute S values - S = Dict{Int,Float64}() + S = Dict{Int64,Float64}() S[9] = C[20]*t^2 + C[21]*t + C[22] S[10] = C[23]*t^2 + C[24]*t + C[25] S[11] = C[26]*t^2 + C[27]*t + C[28] @@ -253,11 +269,7 @@ function calculate_cl(panel::Panel, alpha::Float64) elseif panel.aero_model == "inviscid" return 2π * alpha elseif panel.aero_model == "polar_data" - return linear_interpolation( - panel.polar_data[:,1], - panel.polar_data[:,2]; - extrapolation_bc=Line() - )(alpha) + return panel.polar_data[1](alpha) elseif panel.aero_model == "interpolations" return panel.polar_data[1](alpha, 0.0) else @@ -265,6 +277,10 @@ function calculate_cl(panel::Panel, alpha::Float64) end end +function calculate_cl(polar_data::Nothing, model_type::String, alpha::Float64) + return 2π * alpha +end + """ calculate_cd_cm(panel::Panel, alpha::Float64) @@ -281,19 +297,13 @@ function calculate_cd_cm(panel::Panel, alpha::Float64) elseif panel.aero_model == "inviscid" return 0.0, 0.0 elseif panel.aero_model == "polar_data" - cd = linear_interpolation( - panel.polar_data[:,1], - panel.polar_data[:,3]; - extrapolation_bc=Line() - )(alpha) - cm = linear_interpolation( - panel.polar_data[:,1], - panel.polar_data[:,4]; - extrapolation_bc=Line() - )(alpha) + cd = panel.polar_data[2](alpha) + cm = panel.polar_data[3](alpha) return cd, cm elseif panel.aero_model == "interpolations" - return panel.polar_data[2](alpha, 0.0), panel.polar_data[3](alpha, 0.0) + cd = panel.polar_data[2](alpha, 0.0) + cm = panel.polar_data[3](alpha, 0.0) + return cd, cm else throw(ArgumentError("Unsupported aero model: $(panel.aero_model)")) end diff --git a/src/solver.jl b/src/solver.jl index f76a8399..62031ec1 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -204,17 +204,21 @@ function gamma_loop( v_normal_array = zeros(n_panels) v_tangential_array = zeros(n_panels) - AIC_x, AIC_y, AIC_z = @views wing_aero.AIC[1, :, :], wing_aero.AIC[2, :, :], wing_aero.AIC[3, :, :] - + AIC_x, AIC_y, AIC_z = wing_aero.AIC[1, :, :], wing_aero.AIC[2, :, :], wing_aero.AIC[3, :, :] + + velocity_view_x = @view induced_velocity_all[:, 1] + velocity_view_y = @view induced_velocity_all[:, 2] + velocity_view_z = @view induced_velocity_all[:, 3] + iters = 0 for i in 1:solver.max_iterations iters += 1 gamma .= gamma_new # Calculate induced velocities - mul!(view(induced_velocity_all, :, 1), AIC_x, gamma) - mul!(view(induced_velocity_all, :, 2), AIC_y, gamma) - mul!(view(induced_velocity_all, :, 3), AIC_z, gamma) + mul!(velocity_view_x, AIC_x, gamma) + mul!(velocity_view_y, AIC_y, gamma) + mul!(velocity_view_z, AIC_z, gamma) relative_velocity_array .= va_array .+ induced_velocity_all for i in 1:n_panels @@ -240,7 +244,7 @@ function gamma_loop( end for (i, (panel, alpha)) in enumerate(zip(panels, alpha_array)) - cl_array[i] = calculate_cl(panel, alpha) + cl_array[i] = calculate_cl(panel.polar_data, panel.aero_model, alpha) end gamma_new .= 0.5 .* v_a_array.^2 ./ Umagw_array .* cl_array .* chord_array diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index a0a8d6ff..5c19de70 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -239,9 +239,21 @@ function calculate_new_aero_input(aero_input, 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 = (α, β) -> left_weight * cl_left[α, β] + right_weight * cl_right[α, β] - CD_interp = (α, β) -> left_weight * cd_left[α, β] + right_weight * cd_right[α, β] - CM_interp = (α, β) -> left_weight * cm_left[α, β] + right_weight * cm_right[α, β] + @info "Interpolation object comparison" begin + objects_equal = ( + cl_left === cl_right, + cd_left === cd_right, + cm_left === cm_right + ) + (; objects_equal) + end + + 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[α, β] return ("interpolations", (CL_interp, CD_interp, CM_interp)) From 08c47f3d69b8d2d5129534c19bdf919765aad36d Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sun, 23 Feb 2025 10:40:09 -0600 Subject: [PATCH 13/24] more efficient polars --- examples/ram_air_kite.jl | 72 +++++++++++++++++----------------- src/body_aerodynamics.jl | 8 ++-- src/panel.jl | 65 +++++++++++++++--------------- src/solver.jl | 4 +- src/wing_geometry.jl | 17 ++------ test/test_body_aerodynamics.jl | 4 +- test/test_panel.jl | 11 ------ 7 files changed, 81 insertions(+), 100 deletions(-) diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl index a194a373..4c2a2580 100644 --- a/examples/ram_air_kite.jl +++ b/examples/ram_air_kite.jl @@ -35,17 +35,17 @@ vel_app = [ ] * v_a body_aero.va = vel_app -# Plotting geometry -plot_geometry( - body_aero, - ""; - data_type=".svg", - save_path="", - is_save=false, - is_show=true, - view_elevation=15, - view_azimuth=-120 -) +# # Plotting geometry +# plot_geometry( +# body_aero, +# ""; +# data_type=".svg", +# save_path="", +# is_save=false, +# is_show=true, +# view_elevation=15, +# view_azimuth=-120 +# ) # Solving and plotting distributions @time results = solve(VSM, body_aero) @@ -53,30 +53,30 @@ plot_geometry( CAD_y_coordinates = [panel.aerodynamic_center[2] for panel in body_aero.panels] -plot_distribution( - [CAD_y_coordinates], - [results], - ["VSM"]; - title="CAD_spanwise_distributions_alpha_$(round(aoa, digits=1))_beta_$(round(side_slip, digits=1))_yaw_$(round(yaw_rate, digits=1))_v_a_$(round(v_a, digits=1))", - data_type=".pdf", - is_save=false, - is_show=true -) +# plot_distribution( +# [CAD_y_coordinates], +# [results], +# ["VSM"]; +# title="CAD_spanwise_distributions_alpha_$(round(aoa, digits=1))_beta_$(round(side_slip, digits=1))_yaw_$(round(yaw_rate, digits=1))_v_a_$(round(v_a, digits=1))", +# data_type=".pdf", +# is_save=false, +# is_show=true +# ) -plot_polars( - [VSM], - [body_aero], - [ - "VSM from Ram Air Kite OBJ and DAT file", - ]; - angle_range=range(0, 20, length=20), - angle_type="angle_of_attack", - angle_of_attack=0, - side_slip=0, - v_a=10, - title="ram_kite_panels_$(wing.n_panels)_distribution_$(wing.spanwise_panel_distribution)", - data_type=".pdf", - is_save=false, - is_show=true -) +# plot_polars( +# [VSM], +# [body_aero], +# [ +# "VSM from Ram Air Kite OBJ and DAT file", +# ]; +# angle_range=range(0, 20, length=20), +# angle_type="angle_of_attack", +# angle_of_attack=0, +# side_slip=0, +# v_a=10, +# title="ram_kite_panels_$(wing.n_panels)_distribution_$(wing.spanwise_panel_distribution)", +# data_type=".pdf", +# is_save=false, +# is_show=true +# ) nothing \ No newline at end of file diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 5a4a1745..1d1f08aa 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -226,7 +226,7 @@ function calculate_panel_properties(section_list::Vector{Section}, n_panels::Int end """ - calculate_AIC_matrices(body_aero::BodyAerodynamics, model::String, + calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::String, core_radius_fraction::Float64, va_norm_array::Vector{Float64}, va_unit_array::Matrix{Float64}) @@ -236,7 +236,7 @@ Calculate Aerodynamic Influence Coefficient matrices. Returns: Tuple of (AIC_x, AIC_y, AIC_z) matrices """ -function calculate_AIC_matrices(body_aero::BodyAerodynamics, model::String, +function calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::String, core_radius_fraction::Float64, va_norm_array::Vector{Float64}, va_unit_array::Matrix{Float64}) @@ -249,9 +249,9 @@ function calculate_AIC_matrices(body_aero::BodyAerodynamics, model::String, velocity_induced, tempvel, va_unit, U_2D = zeros(MVec3), zeros(MVec3), zeros(MVec3), zeros(MVec3) # Calculate influence coefficients - for icp in 1:length(body_aero.panels) + for icp in eachindex(body_aero.panels) ep = getproperty(body_aero.panels[icp], evaluation_point) - for jring in 1:body_aero.n_panels + for jring in eachindex(body_aero.panels) va_unit .= @views va_unit_array[jring, :] filaments = body_aero.panels[jring].filaments va_norm = va_norm_array[jring] diff --git a/src/panel.jl b/src/panel.jl index 59989b1d..addd8668 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -36,7 +36,9 @@ mutable struct Panel cl_coefficients::Union{Nothing,Vector{Float64}} cd_coefficients::Union{Nothing,Vector{Float64}} cm_coefficients::Union{Nothing,Vector{Float64}} - polar_data::Union{Nothing, Tuple{Extrapolation}} + cl_interp::Function + cd_interp::Function + cm_interp::Function aerodynamic_center::MVec3 control_point::MVec3 bound_point_1::MVec3 @@ -78,10 +80,8 @@ mutable struct Panel aero_model = isa(section_1.aero_input, String) ? section_1.aero_input : section_1.aero_input[1] # Initialize aerodynamic properties - cl_coeffs = nothing - cd_coeffs = nothing - cm_coeffs = nothing - polar_data = nothing + cl_coeffs, cd_coeffs, cm_coeffs = nothing, nothing, nothing + 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) @@ -92,29 +92,21 @@ mutable struct Panel 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 - ) - polar_data = (cl_interp, cd_interp, cm_interp) + 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 = (α) -> cl_interp_(α) + cd_interp = (α) -> cd_interp_(α) + cm_interp = (α) -> 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 = (α, β) -> 0.5cl_left(α, β) + 0.5cl_right(α, β) - cd_interp = (α, β) -> 0.5cd_left(α, β) + 0.5cd_right(α, β) - cm_interp = (α, β) -> 0.5cm_left(α, β) + 0.5cm_right(α, β) - polar_data = (cl_interp, cd_interp, cm_interp) + cl_interp = cl_left === cl_right ? cl_left : + (α, β) -> 0.5cl_left(α, β) + 0.5cl_right(α, β) + cd_interp = cd_left === cd_right ? cd_left : + (α, β) -> 0.5cd_left(α, β) + 0.5cd_right(α, β) + cm_interp = cm_left === cm_right ? cm_left : + (α, β) -> 0.5cm_left(α, β) + 0.5cm_right(α, β) elseif aero_model != "inviscid" throw(ArgumentError("Unsupported aero model: $aero_model")) end @@ -132,7 +124,8 @@ mutable struct Panel new( TE_point_1, LE_point_1, TE_point_2, LE_point_2, chord, nothing, corner_points, aero_model, - cl_coeffs, cd_coeffs, cm_coeffs, polar_data, + cl_coeffs, cd_coeffs, cm_coeffs, + cl_interp, cd_interp, cm_interp, aerodynamic_center, control_point, bound_point_1, bound_point_2, x_airf, y_airf, z_airf, @@ -269,9 +262,9 @@ function calculate_cl(panel::Panel, alpha::Float64) elseif panel.aero_model == "inviscid" return 2π * alpha elseif panel.aero_model == "polar_data" - return panel.polar_data[1](alpha) + return panel.cl_interp(alpha) elseif panel.aero_model == "interpolations" - return panel.polar_data[1](alpha, 0.0) + return panel.cl_interp(alpha, 0.0) else throw(ArgumentError("Unsupported aero model: $(panel.aero_model)")) end @@ -281,6 +274,14 @@ function calculate_cl(polar_data::Nothing, model_type::String, alpha::Float64) return 2π * alpha end +function calculate_cl(polar_data::Tuple{Function, Function, Function}, model_type::String, alpha::Float64) + if model_type == "polar_data" + return polar_data[1](alpha) + elseif model_type == "interpolations" + return polar_data[1](alpha, 0.0) + end +end + """ calculate_cd_cm(panel::Panel, alpha::Float64) @@ -297,12 +298,12 @@ function calculate_cd_cm(panel::Panel, alpha::Float64) elseif panel.aero_model == "inviscid" return 0.0, 0.0 elseif panel.aero_model == "polar_data" - cd = panel.polar_data[2](alpha) - cm = panel.polar_data[3](alpha) + cd = panel.cd_interp(alpha) + cm = panel.cm_interp(alpha) return cd, cm elseif panel.aero_model == "interpolations" - cd = panel.polar_data[2](alpha, 0.0) - cm = panel.polar_data[3](alpha, 0.0) + cd = panel.cd_interp(alpha, 0.0) + cm = panel.cm_interp(alpha, 0.0) return cd, cm else throw(ArgumentError("Unsupported aero model: $(panel.aero_model)")) diff --git a/src/solver.jl b/src/solver.jl index 95dd9723..d275944f 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -188,7 +188,7 @@ function gamma_loop( relaxation_factor::Float64 ) converged = false - n_panels = body_aero.n_panels + n_panels = length(body_aero.panels) alpha_array = body_aero.alpha_array v_a_array = body_aero.v_a_array Umagw_array = similar(v_a_array) @@ -244,7 +244,7 @@ function gamma_loop( end for (i, (panel, alpha)) in enumerate(zip(panels, alpha_array)) - cl_array[i] = calculate_cl(panel.polar_data, panel.aero_model, alpha) + cl_array[i] = calculate_cl(panel, alpha) end gamma_new .= 0.5 .* v_a_array.^2 ./ Umagw_array .* cl_array .* chord_array diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 04116b73..3a288216 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -239,23 +239,14 @@ function calculate_new_aero_input(aero_input, cl_left, cd_left, cm_left = aero_input[section_index][2] cl_right, cd_right, cm_right = aero_input[section_index + 1][2] - @info "Interpolation object comparison" begin - objects_equal = ( - cl_left === cl_right, - cd_left === cd_right, - cm_left === cm_right - ) - (; objects_equal) - end - - CL_interp = cl_left === cl_right ? cl_left : + cl_interp = cl_left === cl_right ? (α, β) -> cl_left[α, β] : (α, β) -> left_weight * cl_left[α, β] + right_weight * cl_right[α, β] - CD_interp = cd_left === cd_right ? cd_left : + cd_interp = cd_left === cd_right ? (α, β) -> cd_left[α, β] : (α, β) -> left_weight * cd_left[α, β] + right_weight * cd_right[α, β] - CM_interp = cm_left === cm_right ? cm_left : + cm_interp = cm_left === cm_right ? (α, β) -> cm_left[α, β] : (α, β) -> left_weight * cm_left[α, β] + right_weight * cm_right[α, β] - return ("interpolations", (CL_interp, CD_interp, CM_interp)) + return ("interpolations", (cl_interp, cd_interp, cm_interp)) elseif model_type == "lei_airfoil_breukels" tube_diameter_left = aero_input[section_index][2][1] diff --git a/test/test_body_aerodynamics.jl b/test/test_body_aerodynamics.jl index ab6ebc4d..7adf648d 100644 --- a/test/test_body_aerodynamics.jl +++ b/test/test_body_aerodynamics.jl @@ -195,13 +195,13 @@ end va_norm_array = fill(norm(Uinf), length(coord)) va_unit_array = repeat(reshape(Uinf ./ norm(Uinf), 1, 3), length(coord)) calculate_AIC_matrices!( - wing_aero, + body_aero, "VSM", core_radius_fraction, va_norm_array, va_unit_array ) - AIC_x, AIC_y, AIC_z = wing_aero.AIC[1, :, :], wing_aero.AIC[2, :, :], wing_aero.AIC[3, :, :] + AIC_x, AIC_y, AIC_z = body_aero.AIC[1, :, :], body_aero.AIC[2, :, :], body_aero.AIC[3, :, :] # Compare matrices with higher precision for VSM @test isapprox(MatrixU, AIC_x, atol=1e-8) diff --git a/test/test_panel.jl b/test/test_panel.jl index 1085be44..2b7dd5cf 100644 --- a/test/test_panel.jl +++ b/test/test_panel.jl @@ -134,17 +134,6 @@ end # Create panel panel = create_panel(section1, section2) - @testset "Panel Polar Data Initialization" begin - # Test if panel has polar data - @test hasproperty(panel, :polar_data) - @test !isnothing(panel.polar_data) - @test size(panel.polar_data) == size(polar_data) - - # Test if panel polar data is correctly averaged - expected_data = (polar_data + big_polar_data) / 2 - @test isapprox(panel.polar_data, expected_data, rtol=1e-5) - end - @testset "Coefficient Interpolation" begin test_alphas = deg2rad.([-5.0, 0.0, 5.0, 10.0]) for alpha in test_alphas From 96b1b0cc0830e968fe2b59632663efa4c9d08ab8 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sun, 23 Feb 2025 17:43:57 -0600 Subject: [PATCH 14/24] use symbols instead of strings --- examples/code_bench.jl | 138 +++++++++++++++++++++++++++-------------- src/panel.jl | 73 ++++++++++------------ src/solver.jl | 22 +++++-- src/wing_geometry.jl | 44 ++++++------- 4 files changed, 165 insertions(+), 112 deletions(-) diff --git a/examples/code_bench.jl b/examples/code_bench.jl index cfe0f033..f3616a0f 100644 --- a/examples/code_bench.jl +++ b/examples/code_bench.jl @@ -15,24 +15,24 @@ 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") +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") + :inviscid) add_section!(wing, [0.0, -span/2, 0.0], # Right tip LE [chord, -span/2, 0.0], # Right tip TE - "inviscid") + :inviscid) # Step 3: Initialize aerodynamics -wa = WingAerodynamics([wing]) +body_aero = 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 +set_va!(body_aero, (vel_app, 0.0)) # Second parameter is yaw rate # Step 4: Initialize solvers for both LLT and VSM methods vsm_solver = Solver(aerodynamic_model_type="VSM") @@ -40,7 +40,7 @@ vsm_solver = Solver(aerodynamic_model_type="VSM") # Benchmark setup velocity_induced = @MVector zeros(3) # StaticArraysCore.MVector{3, Float64} tempvel = @MVector zeros(3) # StaticArraysCore.MVector{3, Float64} -panel = wa.panels[1] +panel = body_aero.panels[1] ep = @MVector [0.25, 9.5, 0.0] # StaticArraysCore.MVector{3, Float64} evaluation_point_on_bound = true # Bool va_norm = 20.0 # Float64 @@ -51,48 +51,96 @@ gamma = 1.0 # Float64 # Create work vectors tuple of MVector{3, Float64} work_vectors = ntuple(_ -> @MVector(zeros(3)), 10) # NTuple{10, StaticArraysCore.MVector{3, Float64}} -@btime VortexStepMethod.calculate_velocity_induced_single_ring_semiinfinite!( - $velocity_induced, - $tempvel, - $panel.filaments, - $ep, - $evaluation_point_on_bound, - $va_norm, - $va_unit, - $gamma, - $core_radius_fraction, - $work_vectors -) +# @btime VortexStepMethod.calculate_velocity_induced_single_ring_semiinfinite!( +# $velocity_induced, +# $tempvel, +# $panel.filaments, +# $ep, +# $evaluation_point_on_bound, +# $va_norm, +# $va_unit, +# $gamma, +# $core_radius_fraction, +# $work_vectors +# ) -@btime VortexStepMethod.calculate_velocity_induced_single_ring_semiinfinite!( - velocity_induced, - tempvel, - panel.filaments, - ep, - evaluation_point_on_bound, - va_norm, - va_unit, - gamma, - core_radius_fraction, - work_vectors -) +# @btime VortexStepMethod.calculate_velocity_induced_single_ring_semiinfinite!( +# velocity_induced, +# tempvel, +# panel.filaments, +# ep, +# evaluation_point_on_bound, +# va_norm, +# va_unit, +# gamma, +# core_radius_fraction, +# work_vectors +# ) -U_2D = @MVector zeros(3) # StaticArraysCore.MVector{3, Float64} +# U_2D = @MVector zeros(3) # StaticArraysCore.MVector{3, Float64} -@btime VortexStepMethod.calculate_velocity_induced_bound_2D!( - $U_2D, - $panel, - $ep, - $work_vectors -) +# @btime VortexStepMethod.calculate_velocity_induced_bound_2D!( +# $U_2D, +# $panel, +# $ep, +# $work_vectors +# ) + +# model = "VSM" +# n_panels = length(body_aero.panels) +# va_norm_array = zeros(n_panels) +# va_unit_array = zeros(n_panels, 3) +# @btime VortexStepMethod.calculate_AIC_matrices!( +# $body_aero, model, +# $core_radius_fraction, +# $va_norm_array, +# $va_unit_array) -model = "VSM" -va_norm_array = zeros(wa.n_panels) -va_unit_array = zeros(wa.n_panels, 3) -@btime VortexStepMethod.calculate_AIC_matrices!( - $wa, model, - $core_radius_fraction, - $va_norm_array, - $va_unit_array) +n_panels = length(body_aero.panels) +gamma_new = zeros(n_panels) +va_array = zeros(n_panels, 3) +chord_array = zeros(n_panels) +x_airf_array = zeros(n_panels, 3) +y_airf_array = zeros(n_panels, 3) +z_airf_array = zeros(n_panels, 3) +relaxation_factor = 0.5 + +# Fill arrays with data from body_aero +for (i, panel) in enumerate(body_aero.panels) + va_array[i, :] .= panel.va + 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 +end + +# # Benchmark gamma_loop +# @btime VortexStepMethod.gamma_loop( +# $vsm_solver, +# $body_aero, +# $gamma_new, +# $va_array, +# $chord_array, +# $x_airf_array, +# $y_airf_array, +# $z_airf_array, +# $body_aero.panels, +# $relaxation_factor; +# log = false +# ) +# Benchmark gamma_loop +@time VortexStepMethod.gamma_loop( + vsm_solver, + body_aero, + gamma_new, + va_array, + chord_array, + x_airf_array, + y_airf_array, + z_airf_array, + body_aero.panels, + relaxation_factor; + log = false +) nothing \ No newline at end of file diff --git a/src/panel.jl b/src/panel.jl index addd8668..f75288d6 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -11,9 +11,9 @@ Represents a panel in a vortex step method simulation. - `TE_point_2::Vector{MVec3}`: Second trailing edge point - `LE_point_2::Vector{MVec3}`: Second leading edge point - `chord::Float64`: Panel chord length -- `va::Union{Nothing,Vector{Float64}}`: Panel velocity +- `va::Vector{Float64}`: Panel velocity - `corner_points::Matrix{Float64}`: Panel corner points -- `aero_model::String`: Aerodynamic model type +- `aero_model::Symbol`: Aerodynamic model type - `aerodynamic_center::Vector{Float64}`: Panel aerodynamic center - `control_point::Vector{MVec3}`: Panel control point - `bound_point_1::Vector{MVec3}`: First bound point @@ -30,12 +30,12 @@ mutable struct Panel TE_point_2::MVec3 LE_point_2::MVec3 chord::Float64 - va::Union{Nothing,Vector{Float64}} + va::Union{Nothing,MVec3} corner_points::Matrix{Float64} - aero_model::String - cl_coefficients::Union{Nothing,Vector{Float64}} - cd_coefficients::Union{Nothing,Vector{Float64}} - cm_coefficients::Union{Nothing,Vector{Float64}} + aero_model::Symbol + cl_coefficients::Vector{Float64} + cd_coefficients::Vector{Float64} + cm_coefficients::Vector{Float64} cl_interp::Function cd_interp::Function cm_interp::Function @@ -74,18 +74,19 @@ mutable struct Panel corner_points = hcat(LE_point_1, TE_point_1, TE_point_2, LE_point_2) # Validate aero model consistency - if section_1.aero_input[1] != section_2.aero_input[1] + if !(section_1.aero_input === section_2.aero_input) throw(ArgumentError("Both sections must have the same aero_input")) end - aero_model = isa(section_1.aero_input, String) ? section_1.aero_input : section_1.aero_input[1] + aero_model = isa(section_1.aero_input, Symbol) ? section_1.aero_input : section_1.aero_input[1] + # Initialize aerodynamic properties - cl_coeffs, cd_coeffs, cm_coeffs = nothing, nothing, nothing + cl_coeffs, cd_coeffs, cm_coeffs = zeros(3), zeros(3), zeros(3) cl_interp, cd_interp, cm_interp = ()->nothing, ()->nothing, ()->nothing - if aero_model == "lei_airfoil_breukels" + if aero_model === :lei_airfoil_breukels cl_coeffs, cd_coeffs, cm_coeffs = compute_lei_coefficients(section_1, section_2) - elseif aero_model == "polar_data" + 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) @@ -98,16 +99,16 @@ mutable struct Panel cl_interp = (α) -> cl_interp_(α) cd_interp = (α) -> cd_interp_(α) cm_interp = (α) -> cm_interp_(α) - elseif aero_model == "interpolations" + 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 ? cl_left : - (α, β) -> 0.5cl_left(α, β) + 0.5cl_right(α, β) + (α) -> 0.5cl_left(α, 0.0) + 0.5cl_right(α, 0.0) # TODO: add trailing edge deflection cd_interp = cd_left === cd_right ? cd_left : - (α, β) -> 0.5cd_left(α, β) + 0.5cd_right(α, β) + (α) -> 0.5cd_left(α, 0.0) + 0.5cd_right(α, 0.0) cm_interp = cm_left === cm_right ? cm_left : - (α, β) -> 0.5cm_left(α, β) + 0.5cm_right(α, β) - elseif aero_model != "inviscid" + (α) -> 0.5cm_left(α, 0.0) + 0.5cm_right(α, 0.0) + elseif !(aero_model === :inviscid) throw(ArgumentError("Unsupported aero model: $aero_model")) end @@ -252,35 +253,25 @@ Calculate lift coefficient for given angle of attack. # Returns - `Float64`: Lift coefficient (Cl) """ -function calculate_cl(panel::Panel, alpha::Float64) - if panel.aero_model == "lei_airfoil_breukels" +function calculate_cl(panel::Panel, alpha::Float64)::Float64 + cl = 0.0 + if panel.aero_model === :lei_airfoil_breukels cl = evalpoly(rad2deg(alpha), reverse(panel.cl_coefficients)) if abs(alpha) > (π/9) cl = 2 * cos(alpha) * sin(alpha)^2 end - return cl - elseif panel.aero_model == "inviscid" - return 2π * alpha - elseif panel.aero_model == "polar_data" - return panel.cl_interp(alpha) - elseif panel.aero_model == "interpolations" - return panel.cl_interp(alpha, 0.0) + elseif panel.aero_model === :inviscid + cl = 2π * alpha + elseif panel.aero_model === :polar_data + cl = panel.cl_interp(alpha) + elseif panel.aero_model === :interpolations + cl = panel.cl_interp(alpha, 0.0) else throw(ArgumentError("Unsupported aero model: $(panel.aero_model)")) end + return cl end -function calculate_cl(polar_data::Nothing, model_type::String, alpha::Float64) - return 2π * alpha -end - -function calculate_cl(polar_data::Tuple{Function, Function, Function}, model_type::String, alpha::Float64) - if model_type == "polar_data" - return polar_data[1](alpha) - elseif model_type == "interpolations" - return polar_data[1](alpha, 0.0) - end -end """ calculate_cd_cm(panel::Panel, alpha::Float64) @@ -288,20 +279,20 @@ end Calculate drag and moment coefficients for given angle of attack. """ function calculate_cd_cm(panel::Panel, alpha::Float64) - if panel.aero_model == "lei_airfoil_breukels" + if panel.aero_model === :lei_airfoil_breukels cd = evalpoly(rad2deg(alpha), reverse(panel.cd_coefficients)) cm = evalpoly(rad2deg(alpha), reverse(panel.cm_coefficients)) if abs(alpha) > (π/9) # Outside ±20 degrees cd = 2 * sin(alpha)^3 end return cd, cm - elseif panel.aero_model == "inviscid" + elseif panel.aero_model === :inviscid return 0.0, 0.0 - elseif panel.aero_model == "polar_data" + elseif panel.aero_model === :polar_data cd = panel.cd_interp(alpha) cm = panel.cm_interp(alpha) return cd, cm - elseif panel.aero_model == "interpolations" + elseif panel.aero_model === :interpolations cd = panel.cd_interp(alpha, 0.0) cm = panel.cm_interp(alpha, 0.0) return cd, cm diff --git a/src/solver.jl b/src/solver.jl index d275944f..a360954b 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -185,7 +185,8 @@ function gamma_loop( y_airf_array::Matrix{Float64}, z_airf_array::Matrix{Float64}, panels::Vector{Panel}, - relaxation_factor::Float64 + relaxation_factor::Float64; + log::Bool = true ) converged = false n_panels = length(body_aero.panels) @@ -244,7 +245,20 @@ function gamma_loop( end for (i, (panel, alpha)) in enumerate(zip(panels, alpha_array)) - cl_array[i] = calculate_cl(panel, alpha) + @time cl = if panel.aero_model === :lei_airfoil_breukels + cl = evalpoly(rad2deg(alpha), reverse(panel.cl_coefficients)) + if abs(alpha) > (π/9) + cl = 2 * cos(alpha) * sin(alpha)^2 + end + cl + elseif panel.aero_model === :inviscid + 2π * alpha + elseif panel.aero_model === :polar_data + panel.cl_interp(alpha) + elseif panel.aero_model === :interpolations + panel.cl_interp(alpha, 0.0) + end + @time cl_array[i] = calculate_cl(panel, alpha) end gamma_new .= 0.5 .* v_a_array.^2 ./ Umagw_array .* cl_array .* chord_array @@ -276,9 +290,9 @@ function gamma_loop( end end - if converged + if log && converged @info "Converged after $iters iterations" - else + elseif log @warn "NO convergence after $(solver.max_iterations) iterations" end diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 3a288216..6d7dbd60 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -33,7 +33,7 @@ Represents a wing composed of multiple sections with aerodynamic properties. # Fields - `n_panels::Int64`: Number of panels in aerodynamic mesh -- `spanwise_panel_distribution::String`: Panel distribution type +- `spanwise_panel_distribution::Symbol`: Panel distribution type - `spanwise_direction::Vector{Float64}`: Wing span direction vector - `sections::Vector{Section}`: List of wing sections @@ -46,12 +46,12 @@ Represents a wing composed of multiple sections with aerodynamic properties. """ mutable struct Wing <: AbstractWing n_panels::Int64 - spanwise_panel_distribution::String + spanwise_panel_distribution::Symbol spanwise_direction::PosVector sections::Vector{Section} function Wing(n_panels::Int; - spanwise_panel_distribution::String="linear", + spanwise_panel_distribution::Symbol=:linear, spanwise_direction::PosVector=MVec3([0.0, 1.0, 0.0])) new(n_panels, spanwise_panel_distribution, @@ -125,7 +125,7 @@ function refine_aerodynamic_mesh(wing::AbstractWing) end # Handle special cases - if wing.spanwise_panel_distribution == "unchanged" || length(wing.sections) == n_sections + if wing.spanwise_panel_distribution === :unchanged || length(wing.sections) == n_sections return wing.sections end @@ -138,9 +138,9 @@ function refine_aerodynamic_mesh(wing::AbstractWing) end # Handle different distribution types - if wing.spanwise_panel_distribution == "split_provided" + if wing.spanwise_panel_distribution === :split_provided return refine_mesh_by_splitting_provided_sections(wing) - elseif wing.spanwise_panel_distribution in ["linear", "cosine", "cosine_van_Garrel"] + elseif wing.spanwise_panel_distribution in [:linear, :cosine, :cosine_van_Garrel] return refine_mesh_for_linear_cosine_distribution( wing.spanwise_panel_distribution, n_sections, @@ -194,16 +194,16 @@ function calculate_new_aero_input(aero_input, left_weight::Float64, right_weight::Float64) - if aero_input[section_index][1] != aero_input[section_index + 1][1] + if !(aero_input[section_index] === aero_input[section_index + 1]) throw(ArgumentError("Different aero models over the span are not supported")) end - model_type = isa(aero_input[section_index], String) ? aero_input[section_index] : aero_input[section_index][1] + model_type = isa(aero_input[section_index], Symbol) ? aero_input[section_index] : aero_input[section_index][1] - if model_type == "inviscid" - return ("inviscid") + if model_type === :inviscid + return :inviscid - elseif model_type == "polar_data" + elseif model_type === :polar_data polar_left = aero_input[section_index][2] polar_right = aero_input[section_index + 1][2] @@ -233,9 +233,9 @@ function calculate_new_aero_input(aero_input, 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)) + return (:polar_data, hcat(alpha_common, CD_interp, CL_interp, CM_interp)) - elseif model_type == "interpolations" + 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] @@ -246,9 +246,9 @@ function calculate_new_aero_input(aero_input, cm_interp = cm_left === cm_right ? (α, β) -> cm_left[α, β] : (α, β) -> left_weight * cm_left[α, β] + right_weight * cm_right[α, β] - return ("interpolations", (cl_interp, cd_interp, cm_interp)) + return (:interpolations, (cl_interp, cd_interp, cm_interp)) - elseif model_type == "lei_airfoil_breukels" + 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] tube_diameter_i = tube_diameter_left * left_weight + tube_diameter_right * right_weight @@ -260,7 +260,7 @@ function calculate_new_aero_input(aero_input, @debug "Interpolation weights" left_weight right_weight @debug "Interpolated parameters" tube_diameter_i chamber_height_i - return ("lei_airfoil_breukels", [tube_diameter_i, chamber_height_i]) + return (:lei_airfoil_breukels, [tube_diameter_i, chamber_height_i]) else throw(ArgumentError("Unsupported aero model: $(model_type)")) end @@ -268,7 +268,7 @@ end """ refine_mesh_for_linear_cosine_distribution( - spanwise_panel_distribution::String, + spanwise_panel_distribution::Symbol, n_sections::Int, LE::Matrix{Float64}, TE::Matrix{Float64}, @@ -277,7 +277,7 @@ end Refine wing mesh using linear or cosine spacing. # Arguments -- `spanwise_panel_distribution`: Distribution type ("linear", "cosine", or "cosine_van_Garrel") +- `spanwise_panel_distribution`: Distribution type (:linear, :cosine, or :cosine_van_Garrel) - `n_sections`: Number of sections to generate - `LE`: Matrix of leading edge points - `TE`: Matrix of trailing edge points @@ -287,7 +287,7 @@ Returns: Vector{Section}: List of refined sections """ function refine_mesh_for_linear_cosine_distribution( - spanwise_panel_distribution::String, + spanwise_panel_distribution::Symbol, n_sections::Int, LE, TE, @@ -302,9 +302,9 @@ function refine_mesh_for_linear_cosine_distribution( qc_cum_length = vcat(0, cumsum(qc_lengths)) # 2. Define target lengths - target_lengths = if spanwise_panel_distribution == "linear" + target_lengths = if spanwise_panel_distribution === :linear range(0, qc_total_length, n_sections) - elseif spanwise_panel_distribution in ["cosine", "cosine_van_Garrel"] + elseif spanwise_panel_distribution in [:cosine, :cosine_van_Garrel] theta = range(0, π, n_sections) qc_total_length .* (1 .- cos.(theta)) ./ 2 else @@ -368,7 +368,7 @@ function refine_mesh_for_linear_cosine_distribution( end # Apply van Garrel distribution if requested - if spanwise_panel_distribution == "cosine_van_Garrel" + if spanwise_panel_distribution === :cosine_van_Garrel new_sections = calculate_cosine_van_Garrel(new_sections) end From d51848213913a70b393323adb10ed0b08113c14f Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sun, 23 Feb 2025 18:12:44 -0600 Subject: [PATCH 15/24] switch to Symbol --- examples/rectangular_wing.jl | 6 +-- src/body_aerodynamics.jl | 16 +++---- src/kite_geometry.jl | 4 +- src/panel.jl | 32 ++++++------- src/solver.jl | 25 +++------- src/wing_geometry.jl | 18 +++---- test/test_body_aerodynamics.jl | 24 +++++----- test/test_panel.jl | 16 +++---- test/test_plotting.jl | 14 +++--- test/test_wing_geometry.jl | 86 +++++++++++++++++----------------- 10 files changed, 114 insertions(+), 127 deletions(-) diff --git a/examples/rectangular_wing.jl b/examples/rectangular_wing.jl index 3acf8714..f5657030 100644 --- a/examples/rectangular_wing.jl +++ b/examples/rectangular_wing.jl @@ -14,17 +14,17 @@ 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") +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") + :inviscid) add_section!(wing, [0.0, -span/2, 0.0], # Right tip LE [chord, -span/2, 0.0], # Right tip TE - "inviscid") + :inviscid) # Step 3: Initialize aerodynamics wa = BodyAerodynamics([wing]) diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 1d1f08aa..3580c0e8 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -226,7 +226,7 @@ function calculate_panel_properties(section_list::Vector{Section}, n_panels::Int end """ - calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::String, + calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::Symbol, core_radius_fraction::Float64, va_norm_array::Vector{Float64}, va_unit_array::Matrix{Float64}) @@ -236,14 +236,14 @@ Calculate Aerodynamic Influence Coefficient matrices. Returns: Tuple of (AIC_x, AIC_y, AIC_z) matrices """ -function calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::String, +function calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::Symbol, core_radius_fraction::Float64, va_norm_array::Vector{Float64}, va_unit_array::Matrix{Float64}) - model in ["VSM", "LLT"] || throw(ArgumentError("Model must be VSM or LLT")) + model in [:VSM, :LLT] || throw(ArgumentError("Model must be VSM or LLT")) # Determine evaluation point based on model - evaluation_point = model == "VSM" ? :control_point : :aerodynamic_center - evaluation_point_on_bound = model == "LLT" + evaluation_point = model === :VSM ? :control_point : :aerodynamic_center + evaluation_point_on_bound = model === :LLT # Initialize AIC matrices velocity_induced, tempvel, va_unit, U_2D = zeros(MVec3), zeros(MVec3), zeros(MVec3), zeros(MVec3) @@ -270,7 +270,7 @@ function calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::String, body_aero.AIC[:, icp, jring] .= velocity_induced # Subtract 2D induced velocity for VSM - if icp == jring && model == "VSM" + if icp == jring && model === :VSM calculate_velocity_induced_bound_2D!(U_2D, body_aero.panels[jring], ep, body_aero.work_vectors) body_aero.AIC[:, icp, jring] .-= U_2D end @@ -398,7 +398,7 @@ end """ calculate_results(body_aero::BodyAerodynamics, gamma_new::Vector{Float64}, - density::Float64, aerodynamic_model_type::String, + density::Float64, aerodynamic_model_type::Symbol, core_radius_fraction::Float64, mu::Float64, alpha_array::Vector{Float64}, v_a_array::Vector{Float64}, chord_array::Vector{Float64}, x_airf_array::Matrix{Float64}, @@ -415,7 +415,7 @@ Returns: function calculate_results(body_aero::BodyAerodynamics, gamma_new::Vector{Float64}, density::Float64, - aerodynamic_model_type::String, + aerodynamic_model_type::Symbol, core_radius_fraction::Float64, mu::Float64, alpha_array::Vector{Float64}, diff --git a/src/kite_geometry.jl b/src/kite_geometry.jl index 4e53097c..7a2b9f4e 100644 --- a/src/kite_geometry.jl +++ b/src/kite_geometry.jl @@ -176,7 +176,7 @@ Represents a curved wing that inherits from Wing with additional geometric prope # Fields - All fields from Wing: - `n_panels::Int`: Number of panels in aerodynamic mesh - - `spanwise_panel_distribution::String`: Panel distribution type + - `spanwise_panel_distribution::Symbol`: Panel distribution type - `spanwise_direction::Vector{Float64}`: Wing span direction vector - `sections::Vector{Section}`: List of wing sections - Additional fields: @@ -189,7 +189,7 @@ Same as Wing """ mutable struct KiteWing <: AbstractWing n_panels::Int64 - spanwise_panel_distribution::String + spanwise_panel_distribution::Symbol spanwise_direction::Vector{Float64} sections::Vector{Section} diff --git a/src/panel.jl b/src/panel.jl index f75288d6..0caf805d 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -74,15 +74,15 @@ mutable struct Panel corner_points = hcat(LE_point_1, TE_point_1, TE_point_2, LE_point_2) # Validate aero model consistency - if !(section_1.aero_input === section_2.aero_input) - throw(ArgumentError("Both sections must have the same aero_input")) - end aero_model = isa(section_1.aero_input, Symbol) ? section_1.aero_input : section_1.aero_input[1] - + aero_model_2 = isa(section_2.aero_input, Symbol) ? section_2.aero_input : section_2.aero_input[1] + if !(aero_model === aero_model_2) + throw(ArgumentError("Both sections must have the same aero_input, not $aero_model and $aero_model_2")) + end # Initialize aerodynamic properties cl_coeffs, cd_coeffs, cm_coeffs = zeros(3), zeros(3), zeros(3) - cl_interp, cd_interp, cm_interp = ()->nothing, ()->nothing, ()->nothing + cl_interp, cd_interp, cm_interp = (α::Float64)->0.0, (α::Float64)->0.0, (α::Float64)->0.0 if aero_model === :lei_airfoil_breukels cl_coeffs, cd_coeffs, cm_coeffs = compute_lei_coefficients(section_1, section_2) @@ -96,18 +96,18 @@ mutable struct Panel 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 = (α) -> cl_interp_(α) - cd_interp = (α) -> cd_interp_(α) - cm_interp = (α) -> cm_interp_(α) + 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 ? cl_left : - (α) -> 0.5cl_left(α, 0.0) + 0.5cl_right(α, 0.0) # TODO: add trailing edge deflection - cd_interp = cd_left === cd_right ? cd_left : - (α) -> 0.5cd_left(α, 0.0) + 0.5cd_right(α, 0.0) - cm_interp = cm_left === cm_right ? cm_left : - (α) -> 0.5cm_left(α, 0.0) + 0.5cm_right(α, 0.0) + 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) elseif !(aero_model === :inviscid) throw(ArgumentError("Unsupported aero model: $aero_model")) end @@ -263,9 +263,9 @@ function calculate_cl(panel::Panel, alpha::Float64)::Float64 elseif panel.aero_model === :inviscid cl = 2π * alpha elseif panel.aero_model === :polar_data - cl = panel.cl_interp(alpha) + cl = panel.cl_interp(alpha)::Float64 elseif panel.aero_model === :interpolations - cl = panel.cl_interp(alpha, 0.0) + cl = panel.cl_interp(alpha)::Float64 else throw(ArgumentError("Unsupported aero model: $(panel.aero_model)")) end diff --git a/src/solver.jl b/src/solver.jl index a360954b..f2538562 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -6,7 +6,7 @@ Main solver structure for the Vortex Step Method. """ struct Solver # General settings - aerodynamic_model_type::String + aerodynamic_model_type::Symbol density::Float64 max_iterations::Int64 allowed_error::Float64 @@ -18,13 +18,13 @@ struct Solver artificial_damping::NamedTuple{(:k2, :k4), Tuple{Float64, Float64}} # Additional settings - type_initial_gamma_distribution::String + type_initial_gamma_distribution::Symbol core_radius_fraction::Float64 mu::Float64 is_only_f_and_gamma_output::Bool function Solver(; - aerodynamic_model_type::String="VSM", + aerodynamic_model_type::Symbol=:VSM, density::Float64=1.225, max_iterations::Int64=1500, allowed_error::Float64=1e-5, @@ -32,7 +32,7 @@ struct Solver relaxation_factor::Float64=0.03, is_with_artificial_damping::Bool=false, artificial_damping::NamedTuple{(:k2, :k4), Tuple{Float64, Float64}}=(k2=0.1, k4=0.0), - type_initial_gamma_distribution::String="elliptic", + type_initial_gamma_distribution::Symbol=:elliptic, core_radius_fraction::Float64=1e-20, mu::Float64=1.81e-5, is_only_f_and_gamma_output::Bool=false @@ -99,7 +99,7 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n # Initialize gamma distribution gamma_initial = if isnothing(gamma_distribution) - if solver.type_initial_gamma_distribution == "elliptic" + if solver.type_initial_gamma_distribution === :elliptic calculate_circulation_distribution_elliptical_wing(body_aero) else zeros(n_panels) @@ -245,20 +245,7 @@ function gamma_loop( end for (i, (panel, alpha)) in enumerate(zip(panels, alpha_array)) - @time cl = if panel.aero_model === :lei_airfoil_breukels - cl = evalpoly(rad2deg(alpha), reverse(panel.cl_coefficients)) - if abs(alpha) > (π/9) - cl = 2 * cos(alpha) * sin(alpha)^2 - end - cl - elseif panel.aero_model === :inviscid - 2π * alpha - elseif panel.aero_model === :polar_data - panel.cl_interp(alpha) - elseif panel.aero_model === :interpolations - panel.cl_interp(alpha, 0.0) - end - @time cl_array[i] = calculate_cl(panel, alpha) + cl_array[i] = calculate_cl(panel, alpha) end gamma_new .= 0.5 .* v_a_array.^2 ./ Umagw_array .* cl_array .* chord_array diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 6d7dbd60..2f51d7b2 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -38,11 +38,11 @@ Represents a wing composed of multiple sections with aerodynamic properties. - `sections::Vector{Section}`: List of wing sections # Distribution types -- "linear": Linear distribution -- "cosine": Cosine distribution -- "`cosine_van_Garrel`": van Garrel cosine distribution -- "split_provided": Split provided sections -- "unchanged": Keep original sections +- :linear: Linear distribution +- :cosine: Cosine distribution +- :`cosine_van_Garrel`: van Garrel cosine distribution +- :split_provided: Split provided sections +- :unchanged: Keep original sections """ mutable struct Wing <: AbstractWing n_panels::Int64 @@ -194,12 +194,12 @@ function calculate_new_aero_input(aero_input, left_weight::Float64, right_weight::Float64) - if !(aero_input[section_index] === aero_input[section_index + 1]) + model_type = isa(aero_input[section_index], Symbol) ? aero_input[section_index] : aero_input[section_index][1] + model_type_2 = isa(aero_input[section_index + 1], Symbol) ? aero_input[section_index + 1] : aero_input[section_index + 1][1] + if !(model_type === model_type_2) throw(ArgumentError("Different aero models over the span are not supported")) end - model_type = isa(aero_input[section_index], Symbol) ? aero_input[section_index] : aero_input[section_index][1] - if model_type === :inviscid return :inviscid @@ -493,7 +493,7 @@ function refine_mesh_by_splitting_provided_sections(wing::AbstractWing) # Generate sections for this pair new_splitted_sections = refine_mesh_for_linear_cosine_distribution( - "linear", + :linear, num_new_sections + 2, # +2 for endpoints LE_pair, TE_pair, diff --git a/test/test_body_aerodynamics.jl b/test/test_body_aerodynamics.jl index 7adf648d..73285f84 100644 --- a/test/test_body_aerodynamics.jl +++ b/test/test_body_aerodynamics.jl @@ -20,21 +20,21 @@ include("utils.jl") AR = span^2 / (π * span * max_chord / 4) aoa = deg2rad(5) Uinf = [cos(aoa), 0.0, sin(aoa)] .* v_a - model = "VSM" + model = :VSM # Setup wing geometry dist = "cos" core_radius_fraction = 1e-20 coord = generate_coordinates_el_wing(max_chord, span, N, dist) coord_left_to_right = flip_created_coord_in_pairs(deepcopy(coord)) - wing = Wing(N; spanwise_panel_distribution="unchanged") + wing = Wing(N; spanwise_panel_distribution=:unchanged) for idx in 1:2:length(coord_left_to_right[:, 1]) @debug "coord_left_to_right[$idx] = $(coord_left_to_right[idx,:])" add_section!( wing, coord_left_to_right[idx,:], coord_left_to_right[idx+1,:], - "inviscid" + :inviscid ) end @@ -125,13 +125,13 @@ end # Create wing geometry core_radius_fraction = 1e-20 coord_left_to_right = flip_created_coord_in_pairs(deepcopy(coord)) - wing = Wing(n_panels; spanwise_panel_distribution="unchanged") + wing = Wing(n_panels; spanwise_panel_distribution=:unchanged) for idx in 1:2:size(coord_left_to_right, 1) add_section!( wing, coord_left_to_right[idx,:], coord_left_to_right[idx+1,:], - "inviscid" + :inviscid ) end @@ -212,7 +212,7 @@ end @testset "Wing Geometry Creation" begin - function create_geometry(; model="VSM", wing_type="rectangular", plotting=false, N=40) + function create_geometry(; model="VSM", wing_type=:rectangular, plotting=false, N=40) max_chord = 1.0 span = 17.0 AR = span^2 / (π * span * max_chord / 4) @@ -221,7 +221,7 @@ end aoa = 5.7106 * π / 180 Uinf = [cos(aoa), 0.0, sin(aoa)] .* v_a - coord = if wing_type == "rectangular" + coord = if wing_type == :rectangular twist = range(-0.5, 0.5, length=N) beta = range(-2, 2, length=N) generate_coordinates_rect_wing( @@ -232,24 +232,24 @@ end N, "lin" ) - elseif wing_type == "curved" + elseif wing_type == :curved generate_coordinates_curved_wing( max_chord, span, π/4, 5, N, "cos" ) - elseif wing_type == "elliptical" + elseif wing_type == :elliptical generate_coordinates_el_wing(max_chord, span, N, "cos") else error("Invalid wing type") end coord_left_to_right = flip_created_coord_in_pairs(deepcopy(coord)) - wing = Wing(N; spanwise_panel_distribution="unchanged") + wing = Wing(N; spanwise_panel_distribution=:unchanged) for i in 1:2:size(coord_left_to_right, 1) add_section!( wing, coord_left_to_right[i,:], coord_left_to_right[i+1,:], - "inviscid" + :inviscid ) end body_aero = BodyAerodynamics([wing]) @@ -260,7 +260,7 @@ end for model in ["VSM", "LLT"] @debug "model: $model" - for wing_type in ["rectangular", "curved", "elliptical"] + for wing_type in [:rectangular, :curved, :elliptical] @debug "wing_type: $wing_type" body_aero, coord, Uinf, model = create_geometry( model=model, wing_type=wing_type diff --git a/test/test_panel.jl b/test/test_panel.jl index 2b7dd5cf..7cbc795c 100644 --- a/test/test_panel.jl +++ b/test/test_panel.jl @@ -50,8 +50,8 @@ end @testset "Panel Tests" begin @testset "Basic Panel Properties" begin - section1 = Section([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], "inviscid") - section2 = Section([0.0, 10.0, 0.0], [1.0, 10.0, 0.0], "inviscid") + section1 = Section([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], :inviscid) + section2 = Section([0.0, 10.0, 0.0], [1.0, 10.0, 0.0], :inviscid) panel = create_panel(section1, section2) # Test panel initialization @@ -80,8 +80,8 @@ end end @testset "Panel Reference Frame" begin - section1 = Section([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], "inviscid") - section2 = Section([0.0, 10.0, 0.0], [1.0, 10.0, 0.0], "inviscid") + section1 = Section([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], :inviscid) + section2 = Section([0.0, 10.0, 0.0], [1.0, 10.0, 0.0], :inviscid) panel = create_panel(section1, section2) # Test reference frame vectors @@ -91,8 +91,8 @@ end end @testset "Velocity Calculations" begin - section1 = Section([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], "inviscid") - section2 = Section([0.0, 10.0, 0.0], [1.0, 10.0, 0.0], "inviscid") + section1 = Section([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], :inviscid) + section2 = Section([0.0, 10.0, 0.0], [1.0, 10.0, 0.0], :inviscid) panel = create_panel(section1, section2) # Test relative velocity calculations @@ -128,8 +128,8 @@ end end # Create two sections with slightly different polar data - section1 = Section([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], ("polar_data", polar_data)) - section2 = Section([0.0, 10.0, 0.0], [1.0, 10.0, 0.0], ("polar_data", big_polar_data)) + section1 = Section([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], (:polar_data, polar_data)) + section2 = Section([0.0, 10.0, 0.0], [1.0, 10.0, 0.0], (:polar_data, big_polar_data)) # Create panel panel = create_panel(section1, section2) diff --git a/test/test_plotting.jl b/test/test_plotting.jl index 2c348962..daad7556 100644 --- a/test/test_plotting.jl +++ b/test/test_plotting.jl @@ -13,17 +13,17 @@ function create_wa() alpha = deg2rad(alpha_deg) # Step 2: Create wing geometry with linear panel distribution - wing = Wing(n_panels, spanwise_panel_distribution="linear") + 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") + :inviscid) add_section!(wing, [0.0, -span/2, 0.0], # Right tip LE [chord, -span/2, 0.0], # Right tip TE - "inviscid") + :inviscid) # Step 3: Initialize aerodynamics wa = BodyAerodynamics([wing]) @@ -63,8 +63,8 @@ plt.ioff() rm("/tmp/Rectangular_wing_geometry_top_view.pdf") # Step 5: Initialize the solvers - vsm_solver = Solver(aerodynamic_model_type="VSM") - llt_solver = Solver(aerodynamic_model_type="LLT") + vsm_solver = Solver(aerodynamic_model_type=:VSM) + llt_solver = Solver(aerodynamic_model_type=:LLT) # Step 6: Solve the VSM and LLT results_vsm = solve(vsm_solver, wa) @@ -76,7 +76,7 @@ plt.ioff() fig = plot_distribution( [y_coordinates, y_coordinates], [results_vsm, results_llt], - ["VSM", "LLT"], + [:VSM, :LLT], title="Spanwise Distributions" ) @test fig isa plt.PyPlot.Figure @@ -87,7 +87,7 @@ plt.ioff() fig = plot_polars( [llt_solver, vsm_solver], [wa, wa], - ["LLT", "VSM"], + [:LLT, :VSM], angle_range=angle_range, angle_type="angle_of_attack", v_a=v_a, diff --git a/test/test_wing_geometry.jl b/test/test_wing_geometry.jl index a723b44f..25fb9977 100644 --- a/test/test_wing_geometry.jl +++ b/test/test_wing_geometry.jl @@ -19,44 +19,44 @@ end @testset "Wing Geometry Tests" begin @testset "Wing initialization" begin - example_wing = Wing(10; spanwise_panel_distribution="linear") + example_wing = Wing(10; spanwise_panel_distribution=:linear) @test example_wing.n_panels == 10 - @test example_wing.spanwise_panel_distribution == "linear" + @test example_wing.spanwise_panel_distribution == :linear @test example_wing.spanwise_direction ≈ [0.0, 1.0, 0.0] @test length(example_wing.sections) == 0 end @testset "Add section" begin example_wing = Wing(10) - add_section!(example_wing, [0.0, 0.0, 0.0], [-1.0, 0.0, 0.0], "inviscid") + add_section!(example_wing, [0.0, 0.0, 0.0], [-1.0, 0.0, 0.0], :inviscid) @test length(example_wing.sections) == 1 section = example_wing.sections[1] @test section.LE_point ≈ [0.0, 0.0, 0.0] @test section.TE_point ≈ [-1.0, 0.0, 0.0] - @test section.aero_input == "inviscid" + @test section.aero_input == :inviscid end @testset "Robustness left to right" begin example_wing = Wing(10) # Test correct order - add_section!(example_wing, [0.0, 1.0, 0.0], [0.0, 1.0, 0.0], "inviscid") - add_section!(example_wing, [0.0, -1.0, 0.0], [0.0, -1.0, 0.0], "inviscid") - add_section!(example_wing, [0.0, -1.5, 0.0], [0.0, -1.5, 0.0], "inviscid") + add_section!(example_wing, [0.0, 1.0, 0.0], [0.0, 1.0, 0.0], :inviscid) + add_section!(example_wing, [0.0, -1.0, 0.0], [0.0, -1.0, 0.0], :inviscid) + add_section!(example_wing, [0.0, -1.5, 0.0], [0.0, -1.5, 0.0], :inviscid) sections = refine_aerodynamic_mesh(example_wing) # Test right to left order example_wing_1 = Wing(10) - add_section!(example_wing_1, [0.0, -1.5, 0.0], [0.0, -1.5, 0.0], "inviscid") - add_section!(example_wing_1, [0.0, -1.0, 0.0], [0.0, -1.0, 0.0], "inviscid") - add_section!(example_wing_1, [0.0, 1.0, 0.0], [0.0, 1.0, 0.0], "inviscid") + add_section!(example_wing_1, [0.0, -1.5, 0.0], [0.0, -1.5, 0.0], :inviscid) + add_section!(example_wing_1, [0.0, -1.0, 0.0], [0.0, -1.0, 0.0], :inviscid) + add_section!(example_wing_1, [0.0, 1.0, 0.0], [0.0, 1.0, 0.0], :inviscid) sections_1 = refine_aerodynamic_mesh(example_wing_1) # Test random order example_wing_2 = Wing(10) - add_section!(example_wing_2, [0.0, 1.0, 0.0], [0.0, 1.0, 0.0], "inviscid") - add_section!(example_wing_2, [0.0, -1.5, 0.0], [0.0, -1.5, 0.0], "inviscid") - add_section!(example_wing_2, [0.0, -1.0, 0.0], [0.0, -1.0, 0.0], "inviscid") + add_section!(example_wing_2, [0.0, 1.0, 0.0], [0.0, 1.0, 0.0], :inviscid) + add_section!(example_wing_2, [0.0, -1.5, 0.0], [0.0, -1.5, 0.0], :inviscid) + add_section!(example_wing_2, [0.0, -1.0, 0.0], [0.0, -1.0, 0.0], :inviscid) sections_2 = refine_aerodynamic_mesh(example_wing_2) for i in eachindex(sections) @@ -72,9 +72,9 @@ end span = 20.0 # Test linear distribution - wing = Wing(n_panels; spanwise_panel_distribution="linear") - add_section!(wing, [0.0, span/2, 0.0], [-1.0, span/2, 0.0], "inviscid") - add_section!(wing, [0.0, -span/2, 0.0], [-1.0, -span/2, 0.0], "inviscid") + wing = Wing(n_panels; spanwise_panel_distribution=:linear) + add_section!(wing, [0.0, span/2, 0.0], [-1.0, span/2, 0.0], :inviscid) + add_section!(wing, [0.0, -span/2, 0.0], [-1.0, -span/2, 0.0], :inviscid) sections = refine_aerodynamic_mesh(wing) @test length(sections) == wing.n_panels + 1 @@ -87,9 +87,9 @@ end end # Test cosine distribution - wing = Wing(n_panels; spanwise_panel_distribution="cosine") - add_section!(wing, [0.0, span/2, 0.0], [-1.0, span/2, 0.0], "inviscid") - add_section!(wing, [0.0, -span/2, 0.0], [-1.0, -span/2, 0.0], "inviscid") + wing = Wing(n_panels; spanwise_panel_distribution=:cosine) + add_section!(wing, [0.0, span/2, 0.0], [-1.0, span/2, 0.0], :inviscid) + add_section!(wing, [0.0, -span/2, 0.0], [-1.0, -span/2, 0.0], :inviscid) sections = refine_aerodynamic_mesh(wing) @test length(sections) == wing.n_panels + 1 @@ -109,9 +109,9 @@ end n_panels = 1 span = 20.0 - wing = Wing(n_panels; spanwise_panel_distribution="linear") - add_section!(wing, [0.0, span/2, 0.0], [-1.0, span/2, 0.0], "inviscid") - add_section!(wing, [0.0, -span/2, 0.0], [-1.0, -span/2, 0.0], "inviscid") + wing = Wing(n_panels; spanwise_panel_distribution=:linear) + add_section!(wing, [0.0, span/2, 0.0], [-1.0, span/2, 0.0], :inviscid) + add_section!(wing, [0.0, -span/2, 0.0], [-1.0, -span/2, 0.0], :inviscid) sections = refine_aerodynamic_mesh(wing) @test length(sections) == wing.n_panels + 1 @@ -123,9 +123,9 @@ end n_panels = 2 span = 20.0 - wing = Wing(n_panels; spanwise_panel_distribution="linear") - add_section!(wing, [0.0, span/2, 0.0], [-1.0, span/2, 0.0], "inviscid") - add_section!(wing, [0.0, -span/2, 0.0], [-1.0, -span/2, 0.0], "inviscid") + wing = Wing(n_panels; spanwise_panel_distribution=:linear) + add_section!(wing, [0.0, span/2, 0.0], [-1.0, span/2, 0.0], :inviscid) + add_section!(wing, [0.0, -span/2, 0.0], [-1.0, -span/2, 0.0], :inviscid) sections = refine_aerodynamic_mesh(wing) @test length(sections) == wing.n_panels + 1 @@ -138,10 +138,10 @@ end n_panels = 2 span = 20.0 - wing = Wing(n_panels; spanwise_panel_distribution="linear") + wing = Wing(n_panels; spanwise_panel_distribution=:linear) y_coords = [span/2, span/4, 0.0, -span/4, -span/3, -span/2] for y in y_coords - add_section!(wing, [0.0, y, 0.0], [-1.0, y, 0.0], "inviscid") + add_section!(wing, [0.0, y, 0.0], [-1.0, y, 0.0], :inviscid) end sections = refine_aerodynamic_mesh(wing) @@ -159,9 +159,9 @@ end n_panels = 2 span = 10.0 # Total span from -5 to 5 - wing = Wing(n_panels; spanwise_panel_distribution="linear") - add_section!(wing, [0.0, 5.0, 0.0], [-1.0, 5.0, 0.0], "inviscid") - add_section!(wing, [0.0, -5.0, 0.0], [-1.0, -5.0, 0.0], "inviscid") + wing = Wing(n_panels; spanwise_panel_distribution=:linear) + add_section!(wing, [0.0, 5.0, 0.0], [-1.0, 5.0, 0.0], :inviscid) + add_section!(wing, [0.0, -5.0, 0.0], [-1.0, -5.0, 0.0], :inviscid) sections = refine_aerodynamic_mesh(wing) @@ -215,10 +215,10 @@ end n_panels = 4 span = 20.0 - wing = Wing(n_panels; spanwise_panel_distribution="linear") - add_section!(wing, [0.0, span/2, 0.0], [-1.0, span/2, 0.0], ("lei_airfoil_breukels", [0.0, 0.0])) - add_section!(wing, [0.0, 0.0, 0.0], [-1.0, 0.0, 0.0], ("lei_airfoil_breukels", [2.0, 0.5])) - add_section!(wing, [0.0, -span/2, 0.0], [-1.0, -span/2, 0.0], ("lei_airfoil_breukels", [4.0, 1.0])) + wing = Wing(n_panels; spanwise_panel_distribution=:linear) + add_section!(wing, [0.0, span/2, 0.0], [-1.0, span/2, 0.0], (:lei_airfoil_breukels, [0.0, 0.0])) + add_section!(wing, [0.0, 0.0, 0.0], [-1.0, 0.0, 0.0], (:lei_airfoil_breukels, [2.0, 0.5])) + add_section!(wing, [0.0, -span/2, 0.0], [-1.0, -span/2, 0.0], (:lei_airfoil_breukels, [4.0, 1.0])) sections = refine_aerodynamic_mesh(wing) @test length(sections) == wing.n_panels + 1 @@ -234,21 +234,21 @@ end @test isapprox(section.TE_point, expected_TE; rtol=1e-4) aero_input = section.aero_input - @test aero_input[1] == "lei_airfoil_breukels" + @test aero_input[1] == :lei_airfoil_breukels @test isapprox(aero_input[2][1], expected_tube_diameter[i]) @test isapprox(aero_input[2][2], expected_chamber_height[i]) end end @testset "Split provided sections" begin - section1 = Section([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], "inviscid") - section2 = Section([0.0, 1.0, 0.0], [1.0, 1.0, 0.0], "inviscid") - section3 = Section([0.0, 2.0, 0.0], [1.0, 2.0, 0.0], "inviscid") + section1 = Section([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], :inviscid) + section2 = Section([0.0, 1.0, 0.0], [1.0, 1.0, 0.0], :inviscid) + section3 = Section([0.0, 2.0, 0.0], [1.0, 2.0, 0.0], :inviscid) - wing = Wing(6; spanwise_panel_distribution="split_provided") - add_section!(wing, [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], "inviscid") - add_section!(wing, [0.0, 1.0, 0.0], [1.0, 1.0, 0.0], "inviscid") - add_section!(wing, [0.0, 2.0, 0.0], [1.0, 2.0, 0.0], "inviscid") + wing = Wing(6; spanwise_panel_distribution=:split_provided) + add_section!(wing, [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], :inviscid) + add_section!(wing, [0.0, 1.0, 0.0], [1.0, 1.0, 0.0], :inviscid) + add_section!(wing, [0.0, 2.0, 0.0], [1.0, 2.0, 0.0], :inviscid) new_sections = refine_mesh_by_splitting_provided_sections(wing) @@ -262,7 +262,7 @@ end @test 1.0 < new_sections[5].LE_point[2] < 2.0 for section in new_sections - @test section.aero_input == "inviscid" + @test section.aero_input == :inviscid end end end \ No newline at end of file From 7f1d958acbe53d36aafd50406adc79d5ac639d30 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sun, 23 Feb 2025 18:24:03 -0600 Subject: [PATCH 16/24] tests pass again --- src/body_aerodynamics.jl | 19 +++---------------- src/kite_geometry.jl | 4 ++-- test/test_body_aerodynamics.jl | 26 +++++++++++++------------- test/test_kite_geometry.jl | 2 +- test/test_plotting.jl | 4 ++-- test/test_wing_geometry.jl | 10 +++++----- test/thesis_oriol_cayon.jl | 10 +++++----- 7 files changed, 31 insertions(+), 44 deletions(-) diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 3580c0e8..1265c13d 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -375,7 +375,7 @@ function update_effective_angle_of_attack_if_VSM(body_aero::BodyAerodynamics, # Calculate AIC matrices at aerodynamic center using LLT method calculate_AIC_matrices!( - body_aero, "LLT", core_radius_fraction, va_norm_array, va_unit_array + body_aero, :LLT, core_radius_fraction, va_norm_array, va_unit_array ) AIC_x, AIC_y, AIC_z = @views body_aero.AIC[1, :, :], body_aero.AIC[2, :, :], body_aero.AIC[3, :, :] @@ -450,7 +450,7 @@ function calculate_results(body_aero::BodyAerodynamics, moment = reshape((cm_array .* 0.5 .* density .* v_a_array.^2 .* chord_array), :, 1) # Calculate alpha corrections based on model type - alpha_corrected = if aerodynamic_model_type == "VSM" + alpha_corrected = if aerodynamic_model_type === :VSM update_effective_angle_of_attack_if_VSM( body_aero, gamma_new, @@ -461,7 +461,7 @@ function calculate_results(body_aero::BodyAerodynamics, va_norm_array, va_unit_array ) - elseif aerodynamic_model_type == "LLT" + elseif aerodynamic_model_type === :LLT alpha_array else throw(ArgumentError("Unknown aerodynamic model type, should be LLT or VSM")) @@ -541,23 +541,10 @@ function calculate_results(body_aero::BodyAerodynamics, drag_wing_3D_sum += drag_prescribed_va * panel.width side_wing_3D_sum += side_prescribed_va * panel.width - # TODO make this work - # fx_global_3D_sum += fx_global_3D - # fy_global_3D_sum += fy_global_3D - # fz_global_3D_sum += fz_global_3D - # Store coefficients push!(cl_prescribed_va, lift_prescribed_va / (q_inf * panel.chord)) push!(cd_prescribed_va, drag_prescribed_va / (q_inf * panel.chord)) push!(cs_prescribed_va, side_prescribed_va / (q_inf * panel.chord)) - - # TODO translate this - # fx_global_3D_list.append(fx_global_3D) - # fy_global_3D_list.append(fy_global_3D) - # fz_global_3D_list.append(fz_global_3D) - # f_global_3D_list.append( - # np.array([fx_global_3D, fy_global_3D, fz_global_3D]) - # ) end if is_only_f_and_gamma_output diff --git a/src/kite_geometry.jl b/src/kite_geometry.jl index 7a2b9f4e..49070416 100644 --- a/src/kite_geometry.jl +++ b/src/kite_geometry.jl @@ -204,7 +204,7 @@ mutable struct KiteWing <: AbstractWing te_interp::Extrapolation area_interp::Extrapolation - function KiteWing(obj_path, dat_path; alpha=0.0, n_sections=20, crease_frac=0.75, wind_vel=10., mass=1.0, n_panels=54, spanwise_panel_distribution="linear", spanwise_direction=[0.0, 1.0, 0.0]) + function KiteWing(obj_path, dat_path; alpha=0.0, n_sections=20, crease_frac=0.75, wind_vel=10., mass=1.0, n_panels=54, spanwise_panel_distribution=:linear, spanwise_direction=[0.0, 1.0, 0.0]) !isapprox(spanwise_direction, [0.0, 1.0, 0.0]) && @error "Spanwise direction has to be [0.0, 1.0, 0.0]" # Load or create polars @@ -247,7 +247,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 = (:interpolations, (cl_interp, cd_interp, cm_interp)) 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/test/test_body_aerodynamics.jl b/test/test_body_aerodynamics.jl index 73285f84..71ad5924 100644 --- a/test/test_body_aerodynamics.jl +++ b/test/test_body_aerodynamics.jl @@ -140,7 +140,7 @@ end # Calculate reference matrices using thesis functions controlpoints, rings, bladepanels, ringvec, coord_L = - create_geometry_general(coord, Uinf, N, "5fil", "LLT") + create_geometry_general(coord, Uinf, N, "5fil", :LLT) # Test LLT matrices @testset "LLT Matrices" begin @@ -153,7 +153,7 @@ end zeros(N-1), nothing, # data_airf not needed nothing, # conv_crit not needed - "LLT" + :LLT ) # Calculate new matrices @@ -161,7 +161,7 @@ end va_unit_array = repeat(reshape(Uinf ./ norm(Uinf), 1, 3), length(coord)) calculate_AIC_matrices!( body_aero, - "LLT", + :LLT, core_radius_fraction, va_norm_array, va_unit_array @@ -178,7 +178,7 @@ end @testset "VSM Matrices" begin # Calculate reference matrices for VSM controlpoints, rings, bladepanels, ringvec, coord_L = - create_geometry_general(coord, Uinf, N, "5fil", "VSM") + create_geometry_general(coord, Uinf, N, "5fil", :VSM) MatrixU, MatrixV, MatrixW = thesis_induction_matrix_creation( deepcopy(ringvec), @@ -188,7 +188,7 @@ end zeros(N-1), nothing, nothing, - "VSM" + :VSM ) # Calculate new matrices @@ -196,7 +196,7 @@ end va_unit_array = repeat(reshape(Uinf ./ norm(Uinf), 1, 3), length(coord)) calculate_AIC_matrices!( body_aero, - "VSM", + :VSM, core_radius_fraction, va_norm_array, va_unit_array @@ -212,7 +212,7 @@ end @testset "Wing Geometry Creation" begin - function create_geometry(; model="VSM", wing_type=:rectangular, plotting=false, N=40) + function create_geometry(; model=:VSM, wing_type=:rectangular, plotting=false, N=40) max_chord = 1.0 span = 17.0 AR = span^2 / (π * span * max_chord / 4) @@ -221,7 +221,7 @@ end aoa = 5.7106 * π / 180 Uinf = [cos(aoa), 0.0, sin(aoa)] .* v_a - coord = if wing_type == :rectangular + coord = if wing_type === :rectangular twist = range(-0.5, 0.5, length=N) beta = range(-2, 2, length=N) generate_coordinates_rect_wing( @@ -232,11 +232,11 @@ end N, "lin" ) - elseif wing_type == :curved + elseif wing_type === :curved generate_coordinates_curved_wing( max_chord, span, π/4, 5, N, "cos" ) - elseif wing_type == :elliptical + elseif wing_type === :elliptical generate_coordinates_el_wing(max_chord, span, N, "cos") else error("Invalid wing type") @@ -258,7 +258,7 @@ end return body_aero, coord, Uinf, model end - for model in ["VSM", "LLT"] + for model in [:VSM, :LLT] @debug "model: $model" for wing_type in [:rectangular, :curved, :elliptical] @debug "wing_type: $wing_type" @@ -278,7 +278,7 @@ end index_reversed = length(body_aero.panels) - i + 1 panel = body_aero.panels[index_reversed] - evaluation_point = if model == "VSM" + evaluation_point = if model === :VSM panel.control_point else # LLT panel.aerodynamic_center @@ -294,7 +294,7 @@ end atol=1e-4 ) - if model == "VSM" + if model === :VSM @test isapprox( panel.aerodynamic_center, expected_controlpoints[i]["coordinates_aoa"], diff --git a/test/test_kite_geometry.jl b/test/test_kite_geometry.jl index 393809b4..bf962c2c 100644 --- a/test/test_kite_geometry.jl +++ b/test/test_kite_geometry.jl @@ -146,7 +146,7 @@ using Serialization wing = KiteWing(test_obj_path, test_dat_path) @test wing.n_panels == 54 # Default value - @test wing.spanwise_panel_distribution == "linear" + @test wing.spanwise_panel_distribution === :linear @test wing.spanwise_direction ≈ [0.0, 1.0, 0.0] @test length(wing.sections) > 0 # Should have sections now @test wing.mass ≈ 1.0 diff --git a/test/test_plotting.jl b/test/test_plotting.jl index daad7556..92ba03d7 100644 --- a/test/test_plotting.jl +++ b/test/test_plotting.jl @@ -76,7 +76,7 @@ plt.ioff() fig = plot_distribution( [y_coordinates, y_coordinates], [results_vsm, results_llt], - [:VSM, :LLT], + ["VSM", "LLT"], title="Spanwise Distributions" ) @test fig isa plt.PyPlot.Figure @@ -87,7 +87,7 @@ plt.ioff() fig = plot_polars( [llt_solver, vsm_solver], [wa, wa], - [:LLT, :VSM], + ["VSM", "LLT"], angle_range=angle_range, angle_type="angle_of_attack", v_a=v_a, diff --git a/test/test_wing_geometry.jl b/test/test_wing_geometry.jl index 25fb9977..abcaf6b4 100644 --- a/test/test_wing_geometry.jl +++ b/test/test_wing_geometry.jl @@ -14,14 +14,14 @@ Points are compared with approximate equality to handle floating point differenc function ==(a::Section, b::Section) return (isapprox(a.LE_point, b.LE_point; rtol=1e-5, atol=1e-5) && isapprox(a.TE_point, b.TE_point; rtol=1e-5, atol=1e-5) && - a.aero_input == b.aero_input) + a.aero_input === b.aero_input) end @testset "Wing Geometry Tests" begin @testset "Wing initialization" begin example_wing = Wing(10; spanwise_panel_distribution=:linear) @test example_wing.n_panels == 10 - @test example_wing.spanwise_panel_distribution == :linear + @test example_wing.spanwise_panel_distribution === :linear @test example_wing.spanwise_direction ≈ [0.0, 1.0, 0.0] @test length(example_wing.sections) == 0 end @@ -34,7 +34,7 @@ end section = example_wing.sections[1] @test section.LE_point ≈ [0.0, 0.0, 0.0] @test section.TE_point ≈ [-1.0, 0.0, 0.0] - @test section.aero_input == :inviscid + @test section.aero_input === :inviscid end @testset "Robustness left to right" begin @@ -234,7 +234,7 @@ end @test isapprox(section.TE_point, expected_TE; rtol=1e-4) aero_input = section.aero_input - @test aero_input[1] == :lei_airfoil_breukels + @test aero_input[1] === :lei_airfoil_breukels @test isapprox(aero_input[2][1], expected_tube_diameter[i]) @test isapprox(aero_input[2][2], expected_chamber_height[i]) end @@ -262,7 +262,7 @@ end @test 1.0 < new_sections[5].LE_point[2] < 2.0 for section in new_sections - @test section.aero_input == :inviscid + @test section.aero_input === :inviscid end end end \ No newline at end of file diff --git a/test/thesis_oriol_cayon.jl b/test/thesis_oriol_cayon.jl index 23668d6f..f1f096e3 100644 --- a/test/thesis_oriol_cayon.jl +++ b/test/thesis_oriol_cayon.jl @@ -205,7 +205,7 @@ function create_geometry_general(coordinates, Uinf, N, ring_geo, model) normal = x_airf tangential = y_airf - if model == "VSM" + if model === :VSM cp = Dict( "coordinates" => VSMpoint, "chord" => chord, @@ -215,7 +215,7 @@ function create_geometry_general(coordinates, Uinf, N, ring_geo, model) "coordinates_aoa" => LLpoint ) push!(controlpoints, cp) - elseif model == "LLT" + elseif model === :LLT cp = Dict( "coordinates" => LLpoint, "chord" => chord, @@ -320,7 +320,7 @@ function thesis_induction_matrix_creation(ringvec, controlpoints, rings, Uinf, G airf_coord = [controlpoints[icp]["airf_coord"] for icp in 1:N] for icp in 1:N - if model == "VSM" + if model === :VSM # Velocity induced by infinite bound vortex with Gamma = 1 U_2D[icp,:] = velocity_induced_bound_2D(ringvec[icp]) end @@ -347,7 +347,7 @@ function thesis_induction_matrix_creation(ringvec, controlpoints, rings, Uinf, G # Different from old thesis code as this was considered wrong if icp == jring - if model == "VSM" + if model === :VSM MatrixU[icp,jring] -= U_2D[icp,1] MatrixV[icp,jring] -= U_2D[icp,2] MatrixW[icp,jring] -= U_2D[icp,3] @@ -501,7 +501,7 @@ function velocity_induced_single_ring_semiinfinite(ring, controlpoint, model, Ui XV1, Vf, XVP, -GAMMA, Uinf ) elseif filament["id"] == "bound" - if model == "VSM" + if model === :VSM XV2 = filament["x2"] tempvel = velocity_3D_bound_vortex(XV1, XV2, XVP, GAMMA) else From d90015cd56885e022809a3684de396b2d0e310f3 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sun, 23 Feb 2025 19:14:05 -0600 Subject: [PATCH 17/24] working test bench --- examples/code_bench.jl | 126 +++++++++------------ examples/rectangular_wing.jl | 4 +- test/bench.jl | 212 ++++++++++++++++++++++++++++------- 3 files changed, 229 insertions(+), 113 deletions(-) diff --git a/examples/code_bench.jl b/examples/code_bench.jl index f3616a0f..f58b3dca 100644 --- a/examples/code_bench.jl +++ b/examples/code_bench.jl @@ -35,7 +35,7 @@ vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a set_va!(body_aero, (vel_app, 0.0)) # Second parameter is yaw rate # Step 4: Initialize solvers for both LLT and VSM methods -vsm_solver = Solver(aerodynamic_model_type="VSM") +vsm_solver = Solver(aerodynamic_model_type=:VSM) # Benchmark setup velocity_induced = @MVector zeros(3) # StaticArraysCore.MVector{3, Float64} @@ -51,50 +51,50 @@ gamma = 1.0 # Float64 # Create work vectors tuple of MVector{3, Float64} work_vectors = ntuple(_ -> @MVector(zeros(3)), 10) # NTuple{10, StaticArraysCore.MVector{3, Float64}} -# @btime VortexStepMethod.calculate_velocity_induced_single_ring_semiinfinite!( -# $velocity_induced, -# $tempvel, -# $panel.filaments, -# $ep, -# $evaluation_point_on_bound, -# $va_norm, -# $va_unit, -# $gamma, -# $core_radius_fraction, -# $work_vectors -# ) - -# @btime VortexStepMethod.calculate_velocity_induced_single_ring_semiinfinite!( -# velocity_induced, -# tempvel, -# panel.filaments, -# ep, -# evaluation_point_on_bound, -# va_norm, -# va_unit, -# gamma, -# core_radius_fraction, -# work_vectors -# ) - -# U_2D = @MVector zeros(3) # StaticArraysCore.MVector{3, Float64} - -# @btime VortexStepMethod.calculate_velocity_induced_bound_2D!( -# $U_2D, -# $panel, -# $ep, -# $work_vectors -# ) - -# model = "VSM" -# n_panels = length(body_aero.panels) -# va_norm_array = zeros(n_panels) -# va_unit_array = zeros(n_panels, 3) -# @btime VortexStepMethod.calculate_AIC_matrices!( -# $body_aero, model, -# $core_radius_fraction, -# $va_norm_array, -# $va_unit_array) +@btime VortexStepMethod.calculate_velocity_induced_single_ring_semiinfinite!( + $velocity_induced, + $tempvel, + $panel.filaments, + $ep, + $evaluation_point_on_bound, + $va_norm, + $va_unit, + $gamma, + $core_radius_fraction, + $work_vectors +) + +@btime VortexStepMethod.calculate_velocity_induced_single_ring_semiinfinite!( + velocity_induced, + tempvel, + panel.filaments, + ep, + evaluation_point_on_bound, + va_norm, + va_unit, + gamma, + core_radius_fraction, + work_vectors +) + +U_2D = @MVector zeros(3) # StaticArraysCore.MVector{3, Float64} + +@btime VortexStepMethod.calculate_velocity_induced_bound_2D!( + $U_2D, + $panel, + $ep, + $work_vectors +) + +model = :VSM +n_panels = length(body_aero.panels) +va_norm_array = zeros(n_panels) +va_unit_array = zeros(n_panels, 3) +@btime VortexStepMethod.calculate_AIC_matrices!( + $body_aero, model, + $core_radius_fraction, + $va_norm_array, + $va_unit_array) n_panels = length(body_aero.panels) gamma_new = zeros(n_panels) @@ -114,32 +114,18 @@ for (i, panel) in enumerate(body_aero.panels) z_airf_array[i, :] .= panel.z_airf end -# # Benchmark gamma_loop -# @btime VortexStepMethod.gamma_loop( -# $vsm_solver, -# $body_aero, -# $gamma_new, -# $va_array, -# $chord_array, -# $x_airf_array, -# $y_airf_array, -# $z_airf_array, -# $body_aero.panels, -# $relaxation_factor; -# log = false -# ) # Benchmark gamma_loop -@time VortexStepMethod.gamma_loop( - vsm_solver, - body_aero, - gamma_new, - va_array, - chord_array, - x_airf_array, - y_airf_array, - z_airf_array, - body_aero.panels, - relaxation_factor; +@btime VortexStepMethod.gamma_loop( + $vsm_solver, + $body_aero, + $gamma_new, + $va_array, + $chord_array, + $x_airf_array, + $y_airf_array, + $z_airf_array, + $body_aero.panels, + $relaxation_factor; log = false ) diff --git a/examples/rectangular_wing.jl b/examples/rectangular_wing.jl index f5657030..2c432d5f 100644 --- a/examples/rectangular_wing.jl +++ b/examples/rectangular_wing.jl @@ -34,8 +34,8 @@ 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") +llt_solver = Solver(aerodynamic_model_type=:LLT) +vsm_solver = Solver(aerodynamic_model_type=:VSM) # Step 5: Solve using both methods @time results_llt = solve(llt_solver, wa) diff --git a/test/bench.jl b/test/bench.jl index a9a871a6..61f0a3ea 100644 --- a/test/bench.jl +++ b/test/bench.jl @@ -3,8 +3,8 @@ using VortexStepMethod using VortexStepMethod: calculate_AIC_matrices!, gamma_loop, calculate_results, update_effective_angle_of_attack_if_VSM, calculate_projected_area, calculate_cl, calculate_cd_cm, - calculate_velocity_induced_single_ring_semiinfinite, - calculate_velocity_induced_bound_2D, + calculate_velocity_induced_single_ring_semiinfinite!, + calculate_velocity_induced_bound_2D!, velocity_3D_bound_vortex!, velocity_3D_trailing_vortex!, velocity_3D_trailing_vortex_semiinfinite!, @@ -22,22 +22,20 @@ using LinearAlgebra alpha_deg = 30.0 # Angle of attack [degrees] alpha = deg2rad(alpha_deg) - # Create test panels - panels = [] - wing = Wing(n_panels, spanwise_panel_distribution="linear") + 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 - "inviscid") + :inviscid) add_section!(wing, [0.0, -span/2, 0.0], # Right tip LE [chord, -span/2, 0.0], # Right tip TE - "inviscid") + :inviscid) - wing_aero = WingAerodynamics([wing]) + body_aero = BodyAerodynamics([wing]) vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a - set_va!(wing_aero, (vel_app, 0.0)) # Second parameter is yaw rate + set_va!(body_aero, (vel_app, 0.0)) # Second parameter is yaw rate # Initialize solvers for both LLT and VSM methods solver = Solver() @@ -53,65 +51,197 @@ using LinearAlgebra va_norm_array = ones(n_panels) va_unit_array = ones(n_panels, 3) - models = ["VSM", "LLT"] + models = [:VSM, :LLT] core_radius_fractions = [0.001, 10.0] @testset "AIC Matrix Calculation" begin for model in models for frac in core_radius_fractions @testset "Model $model Core Radius Fraction $frac" begin - result = @benchmark calculate_AIC_matrices!($wing_aero, $model, $frac, $va_norm_array, $va_unit_array) - @test result.allocs ≤ 100 # Allow some allocations for matrix setup + global result = @benchmark calculate_AIC_matrices!($body_aero, $model, $frac, $va_norm_array, $va_unit_array) samples = 100 + @test result.allocs ≤ 100 + @info "Model: $(model) \t Core radius fraction: $(frac) \t Allocations: $(result.allocs) \t Memory: $(result.memory)" end end end end - # @testset "Gamma Loop" begin - # result = @benchmark gamma_loop($solver, $wing, $gamma_new, $AIC_x, $AIC_y, $AIC_z) - # @test result.allocs ≤ 50 # Main iteration should be mostly allocation-free - # end + @testset "Gamma Loop" begin + + # Pre-allocate arrays + gamma_new = zeros(n_panels) + va_array = zeros(n_panels, 3) + chord_array = zeros(n_panels) + x_airf_array = zeros(n_panels, 3) + y_airf_array = zeros(n_panels, 3) + z_airf_array = zeros(n_panels, 3) + + # Fill arrays with data + for (i, panel) in enumerate(body_aero.panels) + va_array[i, :] .= panel.va + 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 + end + 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)" + end + end - # @testset "Results Calculation" begin - # result = @benchmark calculate_results($wing, $gamma) - # @test result.allocs ≤ 20 # Allow minimal allocations for results - # end + @testset "Results Calculation" begin + # Pre-allocate arrays + alpha_array = zeros(n_panels) + v_a_array = zeros(n_panels) + chord_array = zeros(n_panels) + x_airf_array = zeros(n_panels, 3) + y_airf_array = zeros(n_panels, 3) + z_airf_array = zeros(n_panels, 3) + va_array = zeros(n_panels, 3) + 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 + + results = @MVector zeros(3) + result = @benchmark calculate_results( + $body_aero, + $gamma, + $density, + :VSM, + 1e-20, + 0.0, + $alpha_array, + $v_a_array, + $chord_array, + $x_airf_array, + $y_airf_array, + $z_airf_array, + $va_array, + $va_norm_array, + $va_unit_array, + $body_aero.panels, + false + ) + @test_broken result.allocs ≤ 100 + end + + # TODO: implement the rest of the benchmarks # @testset "Angle of Attack Update" begin - # result = @benchmark update_effective_angle_of_attack_if_VSM($wing, $gamma) - # @test result.allocs == 0 # Should be allocation-free + # alpha_array = zeros(n_panels) + # result = @benchmark update_effective_angle_of_attack_if_VSM( + # $alpha_array, + # $body_aero, + # $gamma + # ) + # @test result.allocs == 0 # end # @testset "Area Calculations" begin - # result = @benchmark calculate_projected_area($wing) - # @test result.allocs ≤ 10 # Geometric calculations may need some allocations + # area = @MVector zeros(3) + # result = @benchmark calculate_projected_area( + # $area, + # $body_aero + # ) + # @test result.allocs == 0 # end # @testset "Aerodynamic Coefficients" begin - # panel = panels[1] + # panel = body_aero.panels[1] # alpha = 0.1 - # @test (@ballocated calculate_cl($panel, $alpha)) == 0 - # @test (@ballocated calculate_cd_cm($panel, $alpha)) == 0 + # result = @benchmark calculate_cl($panel, $alpha) + # @test result.allocs == 0 + + # cd_cm = @MVector zeros(2) + # result = @benchmark calculate_cd_cm($cd_cm, $panel, $alpha) + # @test result.allocs == 0 # end # @testset "Induced Velocity Calculations" begin + # v_ind = @MVector zeros(3) + # point = @MVector [0.25, 9.5, 0.0] + # work_vectors = ntuple(_ -> @MVector(zeros(3)), 10) + # # Test single ring velocity calculation - # @test (@ballocated calculate_velocity_induced_single_ring_semiinfinite( - # $point, $panels[1], $gamma[1])) == 0 - + # result = @benchmark calculate_velocity_induced_single_ring_semiinfinite!( + # $v_ind, + # $(work_vectors[1]), + # $panel.filaments, + # $point, + # true, + # 20.0, + # $(work_vectors[2]), + # 1.0, + # 1e-20, + # $work_vectors + # ) + # @test result.allocs == 0 + # # Test 2D bound vortex - # @test (@ballocated calculate_velocity_induced_bound_2D( - # $point, $panels[1], $gamma[1])) == 0 - + # result = @benchmark calculate_velocity_induced_bound_2D!( + # $v_ind, + # $panel, + # $point, + # $work_vectors + # ) + # @test result.allocs == 0 + # # Test 3D velocity components - # @test (@ballocated velocity_3D_bound_vortex!( - # $v_ind, $point, $panels[1], $gamma[1])) == 0 - - # @test (@ballocated velocity_3D_trailing_vortex!( - # $v_ind, $point, $panels[1], $gamma[1])) == 0 - - # @test (@ballocated velocity_3D_trailing_vortex_semiinfinite!( - # $v_ind, $point, $panels[1], $gamma[1])) == 0 + # result = @benchmark velocity_3D_bound_vortex!( + # $v_ind, + # $point, + # $panel, + # 1.0, + # 1e-20, + # $work_vectors + # ) + # @test result.allocs == 0 + + # result = @benchmark velocity_3D_trailing_vortex!( + # $v_ind, + # $point, + # $panel, + # 1.0, + # 20.0, + # $work_vectors + # ) + # @test result.allocs == 0 + + # result = @benchmark velocity_3D_trailing_vortex_semiinfinite!( + # $v_ind, + # $point, + # $panel, + # 1.0, + # 20.0, + # $(work_vectors[2]), + # $work_vectors + # ) + # @test result.allocs == 0 # end end \ No newline at end of file From 3ac37398459bd17c774b5fe1ac28f541f92a797b Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sun, 23 Feb 2025 19:24:00 -0600 Subject: [PATCH 18/24] make tests and examples work --- examples/code_bench.jl | 132 ----------------------------------- examples/ram_air_kite.jl | 78 +++++++++++---------- examples/rectangular_wing.jl | 2 +- examples/stall_model.jl | 30 ++++---- src/panel.jl | 22 ++---- test/bench.jl | 7 +- test/runtests.jl | 1 + 7 files changed, 69 insertions(+), 203 deletions(-) delete mode 100644 examples/code_bench.jl diff --git a/examples/code_bench.jl b/examples/code_bench.jl deleted file mode 100644 index f58b3dca..00000000 --- a/examples/code_bench.jl +++ /dev/null @@ -1,132 +0,0 @@ -using Revise -using LinearAlgebra -using ControlPlots -using VortexStepMethod -using StaticArrays -using BenchmarkTools - -# 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 -body_aero = BodyAerodynamics([wing]) - -# Set inflow conditions -vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a -set_va!(body_aero, (vel_app, 0.0)) # Second parameter is yaw rate - -# Step 4: Initialize solvers for both LLT and VSM methods -vsm_solver = Solver(aerodynamic_model_type=:VSM) - -# Benchmark setup -velocity_induced = @MVector zeros(3) # StaticArraysCore.MVector{3, Float64} -tempvel = @MVector zeros(3) # StaticArraysCore.MVector{3, Float64} -panel = body_aero.panels[1] -ep = @MVector [0.25, 9.5, 0.0] # StaticArraysCore.MVector{3, Float64} -evaluation_point_on_bound = true # Bool -va_norm = 20.0 # Float64 -va_unit = @MVector [0.8660254037844387, 0.0, 0.4999999999999999] # StaticArraysCore.MVector{3, Float64} -core_radius_fraction = 1.0e-20 # Float64 -gamma = 1.0 # Float64 - -# Create work vectors tuple of MVector{3, Float64} -work_vectors = ntuple(_ -> @MVector(zeros(3)), 10) # NTuple{10, StaticArraysCore.MVector{3, Float64}} - -@btime VortexStepMethod.calculate_velocity_induced_single_ring_semiinfinite!( - $velocity_induced, - $tempvel, - $panel.filaments, - $ep, - $evaluation_point_on_bound, - $va_norm, - $va_unit, - $gamma, - $core_radius_fraction, - $work_vectors -) - -@btime VortexStepMethod.calculate_velocity_induced_single_ring_semiinfinite!( - velocity_induced, - tempvel, - panel.filaments, - ep, - evaluation_point_on_bound, - va_norm, - va_unit, - gamma, - core_radius_fraction, - work_vectors -) - -U_2D = @MVector zeros(3) # StaticArraysCore.MVector{3, Float64} - -@btime VortexStepMethod.calculate_velocity_induced_bound_2D!( - $U_2D, - $panel, - $ep, - $work_vectors -) - -model = :VSM -n_panels = length(body_aero.panels) -va_norm_array = zeros(n_panels) -va_unit_array = zeros(n_panels, 3) -@btime VortexStepMethod.calculate_AIC_matrices!( - $body_aero, model, - $core_radius_fraction, - $va_norm_array, - $va_unit_array) - -n_panels = length(body_aero.panels) -gamma_new = zeros(n_panels) -va_array = zeros(n_panels, 3) -chord_array = zeros(n_panels) -x_airf_array = zeros(n_panels, 3) -y_airf_array = zeros(n_panels, 3) -z_airf_array = zeros(n_panels, 3) -relaxation_factor = 0.5 - -# Fill arrays with data from body_aero -for (i, panel) in enumerate(body_aero.panels) - va_array[i, :] .= panel.va - 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 -end - -# Benchmark gamma_loop -@btime VortexStepMethod.gamma_loop( - $vsm_solver, - $body_aero, - $gamma_new, - $va_array, - $chord_array, - $x_airf_array, - $y_airf_array, - $z_airf_array, - $body_aero.panels, - $relaxation_factor; - log = false -) - -nothing \ No newline at end of file diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl index 4c2a2580..3accbd39 100644 --- a/examples/ram_air_kite.jl +++ b/examples/ram_air_kite.jl @@ -8,17 +8,19 @@ end using CSV using DataFrames +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", + aerodynamic_model_type=:VSM, is_with_artificial_damping=false ) VSM_with_stall_correction = Solver( - aerodynamic_model_type="VSM", + aerodynamic_model_type=:VSM, is_with_artificial_damping=true ) @@ -35,17 +37,17 @@ vel_app = [ ] * v_a body_aero.va = vel_app -# # Plotting geometry -# plot_geometry( -# body_aero, -# ""; -# data_type=".svg", -# save_path="", -# is_save=false, -# is_show=true, -# view_elevation=15, -# view_azimuth=-120 -# ) +# Plotting geometry +plot && plot_geometry( + body_aero, + ""; + data_type=".svg", + save_path="", + is_save=false, + is_show=true, + view_elevation=15, + view_azimuth=-120 +) # Solving and plotting distributions @time results = solve(VSM, body_aero) @@ -53,30 +55,30 @@ body_aero.va = vel_app CAD_y_coordinates = [panel.aerodynamic_center[2] for panel in body_aero.panels] -# plot_distribution( -# [CAD_y_coordinates], -# [results], -# ["VSM"]; -# title="CAD_spanwise_distributions_alpha_$(round(aoa, digits=1))_beta_$(round(side_slip, digits=1))_yaw_$(round(yaw_rate, digits=1))_v_a_$(round(v_a, digits=1))", -# data_type=".pdf", -# is_save=false, -# is_show=true -# ) +plot && plot_distribution( + [CAD_y_coordinates], + [results], + ["VSM"]; + title="CAD_spanwise_distributions_alpha_$(round(aoa, digits=1))_beta_$(round(side_slip, digits=1))_yaw_$(round(yaw_rate, digits=1))_v_a_$(round(v_a, digits=1))", + data_type=".pdf", + is_save=false, + is_show=true +) -# plot_polars( -# [VSM], -# [body_aero], -# [ -# "VSM from Ram Air Kite OBJ and DAT file", -# ]; -# angle_range=range(0, 20, length=20), -# angle_type="angle_of_attack", -# angle_of_attack=0, -# side_slip=0, -# v_a=10, -# title="ram_kite_panels_$(wing.n_panels)_distribution_$(wing.spanwise_panel_distribution)", -# data_type=".pdf", -# is_save=false, -# is_show=true -# ) +plot && plot_polars( + [VSM], + [body_aero], + [ + "VSM from Ram Air Kite OBJ and DAT file", + ]; + angle_range=range(0, 20, length=20), + angle_type="angle_of_attack", + angle_of_attack=0, + side_slip=0, + v_a=10, + title="ram_kite_panels_$(wing.n_panels)_distribution_$(wing.spanwise_panel_distribution)", + data_type=".pdf", + is_save=false, + is_show=true +) nothing \ No newline at end of file diff --git a/examples/rectangular_wing.jl b/examples/rectangular_wing.jl index 2c432d5f..3ee77fe2 100644 --- a/examples/rectangular_wing.jl +++ b/examples/rectangular_wing.jl @@ -2,7 +2,7 @@ using LinearAlgebra using ControlPlots using VortexStepMethod -plot = false +plot = true # Step 1: Define wing parameters n_panels = 20 # Number of panels diff --git a/examples/stall_model.jl b/examples/stall_model.jl index dd5447fa..49c8cec9 100644 --- a/examples/stall_model.jl +++ b/examples/stall_model.jl @@ -15,7 +15,7 @@ mkpath(save_folder) # Defining discretisation n_panels = 54 -spanwise_panel_distribution = "split_provided" +spanwise_panel_distribution = :split_provided # Load rib data from CSV csv_file_path = joinpath( @@ -30,7 +30,7 @@ rib_list = [] for row in eachrow(df) LE = [row.LE_x, row.LE_y, row.LE_z] TE = [row.TE_x, row.TE_y, row.TE_z] - push!(rib_list, (LE, TE, ("lei_airfoil_breukels", [row.d_tube, row.camber]))) + push!(rib_list, (LE, TE, (:lei_airfoil_breukels, [row.d_tube, row.camber]))) end # Create wing geometry @@ -42,11 +42,11 @@ body_aero_CAD_19ribs = BodyAerodynamics([CAD_wing]) # Create solvers VSM = Solver( - aerodynamic_model_type="VSM", + aerodynamic_model_type=:VSM, is_with_artificial_damping=false ) VSM_with_stall_correction = Solver( - aerodynamic_model_type="VSM", + aerodynamic_model_type=:VSM, is_with_artificial_damping=true ) @@ -63,17 +63,17 @@ vel_app = [ ] * v_a body_aero_CAD_19ribs.va = vel_app -# # Plotting geometry -# plot_geometry( -# body_aero_CAD_19ribs, -# ""; -# data_type=".svg", -# save_path="", -# is_save=false, -# is_show=true, -# view_elevation=15, -# view_azimuth=-120 -# ) +# Plotting geometry +plot_geometry( + body_aero_CAD_19ribs, + ""; + data_type=".svg", + save_path="", + is_save=false, + is_show=true, + view_elevation=15, + view_azimuth=-120 +) # Solving and plotting distributions results = solve(VSM, body_aero_CAD_19ribs) diff --git a/src/panel.jl b/src/panel.jl index 0caf805d..3e813da4 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -262,9 +262,7 @@ function calculate_cl(panel::Panel, alpha::Float64)::Float64 end elseif panel.aero_model === :inviscid cl = 2π * alpha - elseif panel.aero_model === :polar_data - cl = panel.cl_interp(alpha)::Float64 - elseif panel.aero_model === :interpolations + elseif panel.aero_model === :polar_data || panel.aero_model === :interpolations cl = panel.cl_interp(alpha)::Float64 else throw(ArgumentError("Unsupported aero model: $(panel.aero_model)")) @@ -279,26 +277,20 @@ end Calculate drag and moment coefficients for given angle of attack. """ function calculate_cd_cm(panel::Panel, alpha::Float64) + cd, cm = 0.0, 0.0 if panel.aero_model === :lei_airfoil_breukels cd = evalpoly(rad2deg(alpha), reverse(panel.cd_coefficients)) cm = evalpoly(rad2deg(alpha), reverse(panel.cm_coefficients)) if abs(alpha) > (π/9) # Outside ±20 degrees cd = 2 * sin(alpha)^3 end - return cd, cm - elseif panel.aero_model === :inviscid - return 0.0, 0.0 - elseif panel.aero_model === :polar_data - cd = panel.cd_interp(alpha) - cm = panel.cm_interp(alpha) - return cd, cm - elseif panel.aero_model === :interpolations - cd = panel.cd_interp(alpha, 0.0) - cm = panel.cm_interp(alpha, 0.0) - return cd, cm - else + 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 === :inviscid) throw(ArgumentError("Unsupported aero model: $(panel.aero_model)")) end + return cd, cm end """ diff --git a/test/bench.jl b/test/bench.jl index 61f0a3ea..03cbe84c 100644 --- a/test/bench.jl +++ b/test/bench.jl @@ -1,4 +1,5 @@ using BenchmarkTools +using StaticArrays using VortexStepMethod using VortexStepMethod: calculate_AIC_matrices!, gamma_loop, calculate_results, update_effective_angle_of_attack_if_VSM, calculate_projected_area, @@ -55,6 +56,7 @@ using LinearAlgebra core_radius_fractions = [0.001, 10.0] @testset "AIC Matrix Calculation" begin + @info "AIC Matrix Calculation" for model in models for frac in core_radius_fractions @testset "Model $model Core Radius Fraction $frac" begin @@ -67,7 +69,7 @@ using LinearAlgebra end @testset "Gamma Loop" begin - + @info "Gamma Loop" # Pre-allocate arrays gamma_new = zeros(n_panels) va_array = zeros(n_panels, 3) @@ -244,4 +246,5 @@ using LinearAlgebra # ) # @test result.allocs == 0 # end -end \ No newline at end of file +end + diff --git a/test/runtests.jl b/test/runtests.jl index 1a983a64..2864ce96 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Test cd("..") println("Running tests...") @testset verbose = true "Testing VortexStepMethod..." begin + include("bench.jl") include("test_bound_filament.jl") include("test_panel.jl") include("test_semi_infinite_filament.jl") From 49abcf2a883a1ec8e489d6e148a7f7097f811bcd Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Mon, 24 Feb 2025 09:19:16 -0600 Subject: [PATCH 19/24] remove unused dependency --- src/VortexStepMethod.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 7413f616..27483333 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -13,7 +13,6 @@ using NonlinearSolve using Interpolations: linear_interpolation, Line, Extrapolation using Serialization using SharedArrays -using BenchmarkTools # Export public interface export Wing, Section, KiteWing From 1421cfea5d2d84c18650e7c6f01416fd6a8c88f2 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Mon, 24 Feb 2025 16:22:37 +0100 Subject: [PATCH 20/24] add bench.jl --- examples/bench.jl | 46 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 examples/bench.jl diff --git a/examples/bench.jl b/examples/bench.jl new file mode 100644 index 00000000..220f12c3 --- /dev/null +++ b/examples/bench.jl @@ -0,0 +1,46 @@ +using LinearAlgebra +using ControlPlots +using VortexStepMethod + +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) + +# 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 +# time Julia: 0.7 ms + +nothing \ No newline at end of file From 2ebe8969e25dc5ffe41d8fe361cc27ad3d152b17 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Mon, 24 Feb 2025 16:25:42 +0100 Subject: [PATCH 21/24] add bench to the menu --- examples/menu.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/menu.jl b/examples/menu.jl index f7d41d21..c5f635c9 100644 --- a/examples/menu.jl +++ b/examples/menu.jl @@ -3,6 +3,7 @@ using REPL.TerminalMenus options = ["rectangular_wing = include(\"rectangular_wing.jl\")", "ram_air_kite = include(\"ram_air_kite.jl\")", "stall_model = include(\"stall_model.jl\")", + "bench = include(\"bench.jl\")", "quit"] function example_menu() From 9c9d9d939f3a0c5544fb02ccdde8d942688c3496 Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Mon, 24 Feb 2025 16:30:59 +0100 Subject: [PATCH 22/24] less output --- examples/rectangular_wing.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/rectangular_wing.jl b/examples/rectangular_wing.jl index 3ee77fe2..c23a5d51 100644 --- a/examples/rectangular_wing.jl +++ b/examples/rectangular_wing.jl @@ -38,9 +38,9 @@ llt_solver = Solver(aerodynamic_model_type=:LLT) vsm_solver = Solver(aerodynamic_model_type=:VSM) # Step 5: Solve using both methods +results_llt = solve(llt_solver, wa) @time results_llt = solve(llt_solver, wa) -@time results_llt = solve(llt_solver, wa) -@time results_vsm = solve(vsm_solver, wa) +results_vsm = solve(vsm_solver, wa) @time results_vsm = solve(vsm_solver, wa) # time Python: 32.0ms # time Julia: 1.5ms From 60b2a97728dadbc02357ce550bc0ff7bb78f6f4b Mon Sep 17 00:00:00 2001 From: Uwe Fechner Date: Mon, 24 Feb 2025 16:31:55 +0100 Subject: [PATCH 23/24] fix comment --- examples/rectangular_wing.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/rectangular_wing.jl b/examples/rectangular_wing.jl index c23a5d51..61dabc1b 100644 --- a/examples/rectangular_wing.jl +++ b/examples/rectangular_wing.jl @@ -43,7 +43,7 @@ results_llt = solve(llt_solver, wa) results_vsm = solve(vsm_solver, wa) @time results_vsm = solve(vsm_solver, wa) # time Python: 32.0ms -# time Julia: 1.5ms +# time Julia: 0.7ms # Print results comparison println("\nLifting Line Theory Results:") @@ -63,7 +63,6 @@ plot && plot_geometry( is_save=false, is_show=true, ) -nothing # Step 7: Plot spanwise distributions y_coordinates = [panel.aerodynamic_center[2] for panel in wa.panels] From 9be0dd979fb2281ea3c721dce96b107ba855623b Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Mon, 24 Feb 2025 10:08:57 -0600 Subject: [PATCH 24/24] add plot option to example --- examples/stall_model.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/examples/stall_model.jl b/examples/stall_model.jl index 49c8cec9..367513a6 100644 --- a/examples/stall_model.jl +++ b/examples/stall_model.jl @@ -8,6 +8,8 @@ end using CSV using DataFrames +plot = true + # Find root directory root_dir = dirname(@__DIR__) save_folder = joinpath(root_dir, "results", "TUDELFT_V3_LEI_KITE") @@ -64,7 +66,7 @@ vel_app = [ body_aero_CAD_19ribs.va = vel_app # Plotting geometry -plot_geometry( +plot && plot_geometry( body_aero_CAD_19ribs, ""; data_type=".svg", @@ -82,7 +84,7 @@ results = solve(VSM, body_aero_CAD_19ribs) CAD_y_coordinates = [panel.aerodynamic_center[2] for panel in body_aero_CAD_19ribs.panels] -plot_distribution( +plot && plot_distribution( [CAD_y_coordinates, CAD_y_coordinates], [results, results_with_stall], ["VSM", "VSM with stall correction"]; @@ -103,7 +105,7 @@ path_cfd_lebesque = joinpath( "V3_CL_CD_RANS_Lebesque_2024_Rey_300e4.csv" ) -plot_polars( +plot && plot_polars( [VSM, VSM_with_stall_correction], [body_aero_CAD_19ribs, body_aero_CAD_19ribs], [