diff --git a/.gitignore b/.gitignore index 60fe4237..ae204086 100644 --- a/.gitignore +++ b/.gitignore @@ -5,4 +5,3 @@ results/TUDELFT_V3_LEI_KITE/polars/$C_L$ vs $C_D$.pdf *.bin results/TUDELFT_V3_LEI_KITE/polars/tutorial_testing_stall_model_n_panels_54_distribution_split_provided.pdf docs/build/ - 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 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() diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl index a194a373..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 ) @@ -36,7 +38,7 @@ vel_app = [ body_aero.va = vel_app # Plotting geometry -plot_geometry( +plot && plot_geometry( body_aero, ""; data_type=".svg", @@ -53,7 +55,7 @@ plot_geometry( CAD_y_coordinates = [panel.aerodynamic_center[2] for panel in body_aero.panels] -plot_distribution( +plot && plot_distribution( [CAD_y_coordinates], [results], ["VSM"]; @@ -63,7 +65,7 @@ plot_distribution( is_show=true ) -plot_polars( +plot && plot_polars( [VSM], [body_aero], [ diff --git a/examples/rectangular_wing.jl b/examples/rectangular_wing.jl index a73f23ce..61dabc1b 100644 --- a/examples/rectangular_wing.jl +++ b/examples/rectangular_wing.jl @@ -2,6 +2,8 @@ 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] @@ -12,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]) @@ -32,16 +34,16 @@ 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 +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 +# time Julia: 0.7ms # Print results comparison println("\nLifting Line Theory Results:") @@ -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", @@ -61,12 +63,11 @@ 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] -plot_distribution( +plot && plot_distribution( [y_coordinates, y_coordinates], [results_vsm, results_llt], ["VSM", "LLT"], @@ -75,7 +76,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/examples/stall_model.jl b/examples/stall_model.jl index dd5447fa..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") @@ -15,7 +17,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 +32,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 +44,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 +65,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 && 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) @@ -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], [ diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 0182a4ff..bacd950c 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -24,7 +24,8 @@ mutable struct BodyAerodynamics 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 BodyAerodynamics( wings::Vector{T}; @@ -67,8 +68,9 @@ mutable struct BodyAerodynamics 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( panels, wings, @@ -79,7 +81,8 @@ mutable struct BodyAerodynamics stall_angle_list, alpha_array, v_a_array, - work_vectors + work_vectors, + AIC ) end end @@ -223,7 +226,7 @@ function calculate_panel_properties(section_list::Vector{Section}, n_panels::Int end """ - calculate_AIC_matrices(wa::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}) @@ -233,63 +236,65 @@ Calculate Aerodynamic Influence Coefficient matrices. Returns: Tuple of (AIC_x, AIC_y, AIC_z) matrices """ -function calculate_AIC_matrices(wa::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 - AIC = zeros(3, length(wa.panels), length(wa.panels)) + velocity_induced, tempvel, va_unit, U_2D = zeros(MVec3), zeros(MVec3), zeros(MVec3), zeros(MVec3) # Calculate influence coefficients - for icp in 1:length(wa.panels) - ep = getproperty(wa.panels[icp], evaluation_point) - for jring in 1:length(wa.panels) - velocity_induced = calculate_velocity_induced_single_ring_semiinfinite( - wa.panels[jring], + for icp in eachindex(body_aero.panels) + ep = getproperty(body_aero.panels[icp], evaluation_point) + 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] + calculate_velocity_induced_single_ring_semiinfinite!( + velocity_induced, + tempvel, + filaments, ep, evaluation_point_on_bound, - va_norm_array[jring], - va_unit_array[jring, :], + va_norm, + va_unit, 1.0, core_radius_fraction, - wa.work_vectors + body_aero.work_vectors ) - - AIC[:, icp, jring] = velocity_induced + body_aero.AIC[:, icp, jring] .= velocity_induced # Subtract 2D induced velocity for VSM - if icp == jring && model == "VSM" - U_2D = calculate_velocity_induced_bound_2D(wa.panels[jring], ep) - AIC[:, icp, jring] .-= U_2D + 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 end end - - return AIC[1, :, :], AIC[2, :, :], AIC[3, :, :] + return nothing end """ - calculate_circulation_distribution_elliptical_wing(wa::BodyAerodynamics, gamma_0::Float64=1.0) + calculate_circulation_distribution_elliptical_wing(body_aero::BodyAerodynamics, gamma_0::Float64=1.0) Calculate circulation distribution for an elliptical wing. Returns: Vector{Float64}: Circulation distribution along the wing """ -function calculate_circulation_distribution_elliptical_wing(wa::BodyAerodynamics, gamma_0::Float64=1.0) - length(wa.wings) == 1 || throw(ArgumentError("Multiple wings not yet implemented")) +function calculate_circulation_distribution_elliptical_wing(body_aero::BodyAerodynamics, gamma_0::Float64=1.0) + length(body_aero.wings) == 1 || throw(ArgumentError("Multiple wings not yet implemented")) - wing_span = wa.wings[1].span + wing_span = body_aero.wings[1].span @debug "Wing span: $wing_span" # Calculate y-coordinates of control points - y = [panel.control_point[2] for panel in wa.panels] + y = [panel.control_point[2] for panel in body_aero.panels] # Calculate elliptical distribution gamma_i = gamma_0 * sqrt.(1 .- (2 .* y ./ wing_span).^2) @@ -346,7 +351,7 @@ function calculate_stall_angle_list(panels::Vector{Panel}; end """ - update_effective_angle_of_attack_if_VSM(wa::BodyAerodynamics, gamma::Vector{Float64}, + update_effective_angle_of_attack_if_VSM(body_aero::BodyAerodynamics, gamma::Vector{Float64}, core_radius_fraction::Float64, x_airf_array::Matrix{Float64}, y_airf_array::Matrix{Float64}, @@ -359,7 +364,7 @@ Update angle of attack at aerodynamic center for VSM method. Returns: Vector{Float64}: Updated angles of attack """ -function update_effective_angle_of_attack_if_VSM(wa::BodyAerodynamics, +function update_effective_angle_of_attack_if_VSM(body_aero::BodyAerodynamics, gamma::Vector{Float64}, core_radius_fraction::Float64, x_airf_array::Matrix{Float64}, @@ -369,9 +374,10 @@ function update_effective_angle_of_attack_if_VSM(wa::BodyAerodynamics, va_unit_array::Matrix{Float64}) # Calculate AIC matrices at aerodynamic center using LLT method - AIC_x, AIC_y, AIC_z = calculate_AIC_matrices( - wa, "LLT", core_radius_fraction, va_norm_array, va_unit_array + calculate_AIC_matrices!( + 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, :, :] # Calculate induced velocities induced_velocity = [ @@ -391,8 +397,8 @@ function update_effective_angle_of_attack_if_VSM(wa::BodyAerodynamics, end """ - calculate_results(wa::BodyAerodynamics, gamma_new::Vector{Float64}, - density::Float64, aerodynamic_model_type::String, + calculate_results(body_aero::BodyAerodynamics, gamma_new::Vector{Float64}, + 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}, @@ -406,10 +412,10 @@ Calculate final aerodynamic results. Returns: Dict: Results including forces, coefficients and distributions """ -function calculate_results(wa::BodyAerodynamics, +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}, @@ -444,9 +450,9 @@ function calculate_results(wa::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( - wa, + body_aero, gamma_new, core_radius_fraction, x_airf_array, @@ -455,14 +461,14 @@ function calculate_results(wa::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")) end # Verify va is not distributed - if length(wa.va) != 3 + if length(body_aero.va) != 3 throw(ArgumentError("calculate_results not ready for va_distributed input")) end @@ -479,9 +485,9 @@ function calculate_results(wa::BodyAerodynamics, side_wing_3D_sum = 0.0 # Get wing properties - spanwise_direction = wa.wings[1].spanwise_direction - va_mag = norm(wa.va) - va = wa.va + spanwise_direction = body_aero.wings[1].spanwise_direction + va_mag = norm(body_aero.va) + va = body_aero.va va_unit = va / va_mag q_inf = 0.5 * density * va_mag^2 @@ -535,23 +541,10 @@ function calculate_results(wa::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 @@ -562,8 +555,8 @@ function calculate_results(wa::BodyAerodynamics, end # Calculate wing geometry properties - projected_area = sum(wing -> calculate_projected_area(wing), wa.wings) - wing_span = wa.wings[1].span + projected_area = sum(wing -> calculate_projected_area(wing), body_aero.wings) + wing_span = body_aero.wings[1].span aspect_ratio_projected = wing_span^2 / projected_area # Calculate geometric angle of attack @@ -615,16 +608,16 @@ end """ - set_va!(wa::BodyAerodynamics, va::Union{Vector{Float64}, Tuple{Vector{Float64}, Float64}}) + set_va!(body_aero::BodyAerodynamics, va::Union{Vector{Float64}, Tuple{Vector{Float64}, Float64}}) Set velocity array and update wake filaments. # Arguments - `va`: Either a velocity vector or tuple of (velocity vector, yaw_rate) """ -function set_va!(wa::BodyAerodynamics, va) +function set_va!(body_aero::BodyAerodynamics, va) # Add length check for va_vec - if va isa Vector{Float64} && length(va) != 3 && length(va) != length(wa.panels) + if va isa Vector{Float64} && length(va) != 3 && length(va) != length(body_aero.panels) throw(ArgumentError("va must be length 3 or match number of panels")) end # Handle input types @@ -639,15 +632,15 @@ function set_va!(wa::BodyAerodynamics, va) # Calculate va_distribution based on input type va_distribution = if length(va_vec) == 3 && yaw_rate == 0.0 - repeat(reshape(va_vec, 1, 3), length(wa.panels)) - elseif length(va_vec) == length(wa.panels) + repeat(reshape(va_vec, 1, 3), length(body_aero.panels)) + elseif length(va_vec) == length(body_aero.panels) va_vec elseif yaw_rate != 0.0 && length(va_vec) == 3 va_dist = Vector{Float64}[] - for wing in wa.wings + for wing in body_aero.wings # Get spanwise positions - spanwise_positions = [panel.control_point[2] for panel in wa.panels] + spanwise_positions = [panel.control_point[2] for panel in body_aero.panels] # Calculate velocities for each panel for i in 1:wing.n_panels @@ -657,15 +650,15 @@ function set_va!(wa::BodyAerodynamics, va) end reduce(vcat, va_dist) else - throw(ArgumentError("Invalid va distribution: length(va)=$(length(va_vec)) ≠ n_panels=$(length(wa.panels))")) + throw(ArgumentError("Invalid va distribution: length(va)=$(length(va_vec)) ≠ n_panels=$(length(body_aero.panels))")) end # Update panel velocities - for (i, panel) in enumerate(wa.panels) + for (i, panel) in enumerate(body_aero.panels) panel.va = va_distribution[i,:] end # Update wake elements - wa.panels = frozen_wake(va_distribution, wa.panels) - wa._va = va_vec + body_aero.panels = frozen_wake(va_distribution, body_aero.panels) + body_aero._va = va_vec end diff --git a/src/filament.jl b/src/filament.jl index 41730627..35ec2dfa 100644 --- a/src/filament.jl +++ b/src/filament.jl @@ -37,12 +37,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 @@ -61,10 +61,10 @@ 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 - @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) @@ -98,13 +98,13 @@ 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!( - vel::VelVector, +@inline function velocity_3D_trailing_vortex!( + 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 @@ -128,7 +128,7 @@ function velocity_3D_trailing_vortex!( 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) + @@ -166,13 +166,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 @@ -188,7 +188,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/kite_geometry.jl b/src/kite_geometry.jl index 511a56bd..84e088c5 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::MVec3`: 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::MVec3 sections::Vector{Section} @@ -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/src/panel.jl b/src/panel.jl index 78aee0c0..999de0dd 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -13,7 +13,7 @@ Represents a panel in a vortex step method simulation. - `chord::Float64`: Panel chord length - `va::Union{Nothing, MVec3}`: 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 @@ -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 @@ -32,11 +32,13 @@ mutable struct Panel{P} chord::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}} - polar_data::P + aero_model::Symbol + cl_coefficients::Vector{Float64} + cd_coefficients::Vector{Float64} + cm_coefficients::Vector{Float64} + cl_interp::Function + cd_interp::Function + cm_interp::Function aerodynamic_center::MVec3 control_point::MVec3 bound_point_1::MVec3 @@ -72,34 +74,41 @@ mutable struct Panel{P} 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] - throw(ArgumentError("Both sections must have the same aero_input")) + 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 - 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 = zeros(3), zeros(3), zeros(3) + cl_interp, cd_interp, cm_interp = (α::Float64)->0.0, (α::Float64)->0.0, (α::Float64)->0.0 - 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) throw(ArgumentError("Polar data must have same shape")) end polar_data = (aero_1 + aero_2) / 2 - elseif aero_model == "interpolations" + cl_interp_ = linear_interpolation(polar_data[:,1], polar_data[:,2]; extrapolation_bc=NaN) + cd_interp_ = linear_interpolation(polar_data[:,1], polar_data[:,3]; extrapolation_bc=NaN) + cm_interp_ = linear_interpolation(polar_data[:,1], polar_data[:,4]; extrapolation_bc=NaN) + cl_interp = (α::Float64) -> cl_interp_(α) + cd_interp = (α::Float64) -> cd_interp_(α) + cm_interp = (α::Float64) -> cm_interp_(α) + elseif aero_model === :interpolations cl_left, cd_left, cm_left = section_1.aero_input[2] cl_right, cd_right, cm_right = section_2.aero_input[2] - cl_interp = (α, β) -> 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) - elseif aero_model != "inviscid" + 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 @@ -113,10 +122,11 @@ 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, + 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, @@ -243,60 +253,44 @@ 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 linear_interpolation( - panel.polar_data[:,1], - panel.polar_data[:,2]; - extrapolation_bc=Line() - )(alpha) - elseif panel.aero_model == "interpolations" - return panel.polar_data[1](alpha, 0.0) + elseif panel.aero_model === :inviscid + cl = 2π * alpha + 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)")) end + return cl end + """ calculate_cd_cm(panel::Panel, alpha::Float64) 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" + 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 = 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) - return cd, cm - elseif panel.aero_model == "interpolations" - return panel.polar_data[2](alpha, 0.0), panel.polar_data[3](alpha, 0.0) - 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 """ @@ -338,7 +332,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,49 +355,50 @@ 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( - panel::Panel, - evaluation_point::PosVector, +@inline function calculate_velocity_induced_single_ring_semiinfinite!( + velind::MVec3, + tempvel::MVec3, + filaments::Vector{Union{VortexStepMethod.BoundFilament, VortexStepMethod.SemiInfiniteFilament}}, + evaluation_point::MVec3, evaluation_point_on_bound::Bool, va_norm::Float64, - va_unit::VelVector, + va_unit::MVec3, gamma::Float64, core_radius_fraction::Float64, - work_vectors::NTuple{10,Vector{Float64}} + work_vectors::NTuple{10, MVec3} ) - velind = zeros(Float64, 3) - tempvel = zeros(3) - + 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 .= zeros(Float64, 3) + 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, + filaments[i], evaluation_point, gamma, va_norm, work_vectors ) - elseif i ∈ (4, 5) # semi-infinite trailing filaments + elseif i == 4 || i == 5 # semi-infinite trailing filaments velocity_3D_trailing_vortex_semiinfinite!( tempvel, - filament, + filaments[i], va_unit, evaluation_point, gamma, @@ -409,16 +406,15 @@ function calculate_velocity_induced_single_ring_semiinfinite( work_vectors ) else - tempvel .= zeros(Float64, 3) + tempvel .= 0.0 end velind .+= tempvel end - - return velind + return nothing 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 +426,24 @@ 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) + 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/src/solver.jl b/src/solver.jl index 9234c575..5ff23490 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, # TODO do not use magic constants is_only_f_and_gamma_output::Bool = false @@ -74,7 +74,7 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, 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 @@ -89,7 +89,7 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, 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!( body_aero, solver.aerodynamic_model_type, solver.core_radius_fraction, @@ -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) @@ -116,9 +116,6 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n solver, body_aero, gamma_initial, - AIC_x, - AIC_y, - AIC_z, va_array, chord_array, x_airf_array, @@ -134,9 +131,6 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n solver, body_aero, gamma_initial, - AIC_x, - AIC_y, - AIC_z, va_array, chord_array, x_airf_array, @@ -170,6 +164,8 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, 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}, @@ -183,64 +179,94 @@ function gamma_loop( solver::Solver, body_aero::BodyAerodynamics, 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}, 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) alpha_array = body_aero.alpha_array v_a_array = body_aero.v_a_array + Umagw_array = similar(v_a_array) + + gamma = copy(gamma_new) + abs_gamma_new = copy(gamma_new) + 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) + cl_array = zeros(n_panels) + damp = zeros(length(gamma)) + v_normal_array = zeros(n_panels) + v_tangential_array = zeros(n_panels) + + AIC_x, AIC_y, AIC_z = body_aero.AIC[1, :, :], body_aero.AIC[2, :, :], body_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 = copy(gamma_new) + gamma .= gamma_new # Calculate induced velocities - induced_velocity_all = hcat( - AIC_x * gamma, - AIC_y * gamma, - 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 - 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 v_a_array[i] = norm(relative_velocity_crossz[i, :]) + @views Umagw_array[i] = norm(Uinfcrossz_array[i, :]) + end - v_a_array .= norm.(relative_velocity_crossz) - v_aw_array = norm.(Uinfcrossz_array) - - cl_array = [calculate_cl(panel, alpha) for (panel, alpha) in zip(panels, alpha_array)] - gamma_new .= 0.5 .* v_a_array.^2 ./ v_aw_array .* cl_array .* chord_array + for (i, (panel, alpha)) in enumerate(zip(panels, alpha_array)) + cl_array[i] = calculate_cl(panel, alpha) + end + gamma_new .= 0.5 .* v_a_array.^2 ./ Umagw_array .* cl_array .* chord_array # Apply damping if needed if solver.is_with_artificial_damping 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" @@ -251,9 +277,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 4a71deff..45f015f8 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -33,25 +33,25 @@ 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 # 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 - 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] + 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], String) ? 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,19 +233,22 @@ 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] - 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[α, β] + 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)) + 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 @@ -257,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 @@ -265,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}, @@ -274,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 @@ -284,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, @@ -299,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 @@ -365,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 @@ -490,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/bench.jl b/test/bench.jl new file mode 100644 index 00000000..03cbe84c --- /dev/null +++ b/test/bench.jl @@ -0,0 +1,250 @@ +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, + 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) + + 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) + + body_aero = BodyAerodynamics([wing]) + + vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a + set_va!(body_aero, (vel_app, 0.0)) # Second parameter is yaw rate + + # 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 + @info "AIC Matrix Calculation" + for model in models + for frac in core_radius_fractions + @testset "Model $model Core Radius Fraction $frac" begin + global result = @benchmark calculate_AIC_matrices!($body_aero, $model, $frac, $va_norm_array, $va_unit_array) samples = 100 + @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 + @info "Gamma Loop" + # 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 + # 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 + # 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 + # area = @MVector zeros(3) + # result = @benchmark calculate_projected_area( + # $area, + # $body_aero + # ) + # @test result.allocs == 0 + # end + + # @testset "Aerodynamic Coefficients" begin + # panel = body_aero.panels[1] + # alpha = 0.1 + + # 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 + # 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 + # result = @benchmark calculate_velocity_induced_bound_2D!( + # $v_ind, + # $panel, + # $point, + # $work_vectors + # ) + # @test result.allocs == 0 + + # # Test 3D velocity components + # 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 + 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") diff --git a/test/test_body_aerodynamics.jl b/test/test_body_aerodynamics.jl index e2dc2fe8..71ad5924 100644 --- a/test/test_body_aerodynamics.jl +++ b/test/test_body_aerodynamics.jl @@ -1,6 +1,6 @@ # using VortexStepMethod: Wing, BodyAerodynamics, 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 @@ -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 @@ -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,19 +153,20 @@ end zeros(N-1), nothing, # data_airf not needed nothing, # conv_crit not needed - "LLT" + :LLT ) # 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!( body_aero, - "LLT", + :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, :, :] # Compare matrices @test isapprox(MatrixU, AIC_x, atol=1e-5) @@ -177,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), @@ -187,19 +188,20 @@ end zeros(N-1), nothing, nothing, - "VSM" + :VSM ) # 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!( body_aero, - "VSM", + :VSM, core_radius_fraction, va_norm_array, va_unit_array ) + 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) @@ -210,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) @@ -219,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( @@ -230,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]) @@ -256,9 +258,9 @@ 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"] + 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 @@ -276,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 @@ -292,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_panel.jl b/test/test_panel.jl index 1085be44..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,23 +128,12 @@ 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) - @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 diff --git a/test/test_plotting.jl b/test/test_plotting.jl index 2c348962..92ba03d7 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) @@ -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 a723b44f..abcaf6b4 100644 --- a/test/test_wing_geometry.jl +++ b/test/test_wing_geometry.jl @@ -14,49 +14,49 @@ 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") + 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 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