cant remove calculate_velocity_induced_single_ring_semiinfinite! allocs

Author: 1-Bart-1

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

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
 
 """

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,14 +195,16 @@ 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)
     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 = @views wing_aero.AIC[1, :, :], wing_aero.AIC[2, :, :], wing_aero.AIC[3, :, :]
 
     iters = 0
     for i in 1:solver.max_iterations

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 = [

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

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)