Reduce allocation (#126)

* Reduce allocations

* Use normalize!

* Reduce allocations

* Use \cdot

* Cleanup

* Replace cross with \times

* Add bench2.jl

* Update comments

* Add P to VSMSolution and bench2.jl

* Make use of panel_width_array

* Now 112, 600 allocations

* Now 596 allocations

* Down to 590 allocations

* Cleanup

* Update comment

* Fix tests

* Fix examples

* Cleanup

* Down to 586 allocations

* Refactor update_effective_angle_of_attack_if_VSM()

* Now 424 allocations

* Now 422 allocations

* Now 420 allocations

* Change comment

* Cleanup

* Add PreallocationTools

* Now 408 allocations

* Now 398 allocations

* Refactoring

* Refactoring

* Refactoring

* Refactoring

* Cleanup

* Refactoring

* Refactoring

* Update comment

* Now 394 allocations

* Refactoring

* Add BaseResult

* Improve LoopResult

* Refactoring, use LoopResult

* Refactoring

* Add comment

* Refactoring

* Refactoring

* Update docstring

* Refactoring

* Cleanup

* Refactoring

* Refactoring

* Refactoring

* Refactoring

* Update struct BaseResult

* Add minimal working example

* Improve mwe_01.jl

* Now 390 allocations

* Refactoring

* Add TODO

* Now 350 allocations

* Now 348 allocations

* Now 346 allocations

* Now 344 allocations

* Now 342 allocations

* Add comment

* Refactoring

* Refactoring

* Update comment

* Update comments

* Refactoring

* Now 340 allocations

* Add comment

* Now 337 allocations

* Update comments

* Now 328 allocations

* Cleanup

* Update comment

* Update comment

* Update comment

* Add mwe_02.jl

* Cleanup

* Add file to .gitignore

* Use Julia 1.10 to create bin files

* Address remarks

* Remove print statement

* Add comments

---------

Co-authored-by: Uwe Fechner <[email protected]>
Co-authored-by: 1-Bart-1 <[email protected]>

Author: ufechner7

Manifest-v1.10.toml

@@ -15,6 +15,7 @@ Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -47,6 +48,7 @@ Measures = "0.3"
 NonlinearSolve = "4"
 Parameters = "0.12"
 Pkg = "1"
+PreallocationTools = "0.4.25"
 Serialization = "1"
 SharedArrays = "1"
 StaticArrays = "1"

Project.toml

@@ -39,18 +39,18 @@ vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
 set_va!(wa, vel_app)
 
 # Step 4: Initialize solvers for both LLT and VSM methods
-llt_solver = Solver(aerodynamic_model_type=LLT)
-vsm_solver = Solver(aerodynamic_model_type=VSM)
+P = length(wa.panels)
+llt_solver = Solver{P}(aerodynamic_model_type=LLT)
+vsm_solver = Solver{P}(aerodynamic_model_type=VSM)
 
 # Step 5: Solve using both methods
 results_vsm = solve(vsm_solver, wa)
 sol = solve!(vsm_solver, wa)
-results_vsm_base = solve_base(vsm_solver, wa)
-println("Rectangular wing, solve_base:")
-@time results_vsm_base = solve_base(vsm_solver, wa)
-# time Python: 32.0 ms Ryzen 7950x
-# time Julia:   0.6 ms Ryzen 7950x
-#               0.8 ms laptop, performance mode, battery 
+results_vsm_base = solve_base!(vsm_solver, wa)
+println("Rectangular wing, solve_base!:")
+@time results_vsm_base = solve_base!(vsm_solver, wa)
+# time Python: 32.0  ms Ryzen 7950x
+# time Julia:   0.42 ms Ryzen 7950x
 println("Rectangular wing, solve!:")
 @time sol = solve!(vsm_solver, wa)
 println("Rectangular wing, solve:")
@@ -61,7 +61,8 @@ wing = RamAirWing("data/ram_air_kite_body.obj", "data/ram_air_kite_foil.dat")
 body_aero = BodyAerodynamics([wing])
 
 # Create solvers
-vsm_solver = Solver(
+P = length(wa.panels)
+vsm_solver = Solver{P}(
     aerodynamic_model_type=VSM,
     is_with_artificial_damping=false
 )
@@ -81,8 +82,8 @@ set_va!(body_aero, vel_app)
 
 # Solving and plotting distributions
 results = solve(vsm_solver, body_aero)
-results_base = solve_base(vsm_solver, body_aero)
+results_base = solve_base!(vsm_solver, body_aero)
 println("RAM-air kite:")
-@time results_base = solve_base(vsm_solver, body_aero)
+@time results_base = solve_base!(vsm_solver, body_aero)
 
 nothing

examples/bench.jl

@@ -24,7 +24,8 @@ if DEFORM
 end
 
 # Create solvers
-vsm_solver = Solver(
+P = length(body_aero.panels)
+vsm_solver = Solver{P}(
     aerodynamic_model_type=VSM,
     is_with_artificial_damping=false
 )

examples/ram_air_kite.jl

@@ -35,8 +35,9 @@ vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
 set_va!(wa, vel_app, [0, 0, 0.1])
 
 # Step 4: Initialize solvers for both LLT and VSM methods
-llt_solver = Solver(aerodynamic_model_type=LLT)
-vsm_solver = Solver(aerodynamic_model_type=VSM)
+P = length(wa.panels)
+llt_solver = Solver{P}(aerodynamic_model_type=LLT)
+vsm_solver = Solver{P}(aerodynamic_model_type=VSM)
 
 # Step 5: Solve using both methods
 results_llt = solve(llt_solver, wa)

examples/rectangular_wing.jl

@@ -45,11 +45,12 @@ end
 body_aero_CAD_19ribs = BodyAerodynamics([CAD_wing])
 
 # Create solvers
-vsm_solver = Solver(
+P = length(wa.panels)
+vsm_solver = Solver{P}(
     aerodynamic_model_type=VSM,
     is_with_artificial_damping=false
 )
-VSM_with_stall_correction = Solver(
+VSM_with_stall_correction = Solver{P}(
     aerodynamic_model_type=VSM,
     is_with_artificial_damping=true
 )

examples/stall_model.jl

@@ -0,0 +1,63 @@
+# Replace va_norm_array = norm.(eachrow(solver.sol.va_array)) with a for loop
+
+# Testcase that shows that the new function is equivalent to the old, allocating line of code.
+using Test
+using LinearAlgebra  # for the norm function
+
+# Define a mock Solver struct
+struct MockSolver
+    sol::NamedTuple
+end
+
+function calc_norm_array!(va_norm_array, va_array)
+    for i in 1:size(va_array, 1)
+        va_norm_array[i] = norm(view(va_array, i, :))
+    end
+end
+
+@testset "va_norm_array calculation" begin
+    global va_norm_array
+
+    # Create a sample 2D array
+    sample_va_array = [
+        1.0 2.0 3.0;
+        4.0 5.0 6.0;
+        7.0 8.0 9.0
+    ]
+
+    # Create a mock solver with the sample array
+    mock_solver = MockSolver((va_array = sample_va_array,))
+
+    # Calculate va_norm_array
+    n = @allocated va_norm_array = norm.(eachrow(mock_solver.sol.va_array))
+    println(n)
+
+    va_norm_array2 = zeros(3)
+    m = @allocated calc_norm_array!(va_norm_array2, sample_va_array)
+    println(m)
+
+    # Expected results (calculated manually)
+    expected_norms = [
+        sqrt(1^2 + 2^2 + 3^2),
+        sqrt(4^2 + 5^2 + 6^2),
+        sqrt(7^2 + 8^2 + 9^2)
+    ]
+
+    # Test the results
+    @test length(va_norm_array) == size(sample_va_array, 1)
+    @test va_norm_array ≈ expected_norms atol=1e-10
+
+    # Test individual values
+    @test va_norm_array[1] ≈ norm(sample_va_array[1, :]) atol=1e-10
+    @test va_norm_array[2] ≈ norm(sample_va_array[2, :]) atol=1e-10
+    @test va_norm_array[3] ≈ norm(sample_va_array[3, :]) atol=1e-10
+
+    @test length(va_norm_array2) == size(sample_va_array, 1)
+    @test va_norm_array2 ≈ expected_norms atol=1e-10
+
+    # Test individual values
+    @test va_norm_array2[1] ≈ norm(sample_va_array[1, :]) atol=1e-10
+    @test va_norm_array2[2] ≈ norm(sample_va_array[2, :]) atol=1e-10
+    @test va_norm_array2[3] ≈ norm(sample_va_array[3, :]) atol=1e-10
+end
+nothing

mwes/mwe_01.jl

@@ -0,0 +1,14 @@
+
+# I tried to put a lazy cache (to avoid allocations for temporary vectors)
+# into a struct. This does not work at all with the macro provided by the package 
+# Parameters. It works (no syntax error) with the macro Base.@kwdef, but if you 
+# make use of this cache it is extremely slow.
+#
+# This mwe can serve as starting point for further investigations. Currently the only
+# cache that works is a cache in a global const variable.
+
+using PreallocationTools, Parameters
+
+Base.@kwdef mutable struct MockSolver
+    const ca::NTuple{11, LazyBufferCache{typeof(identity), typeof(identity)}} = ([LazyBufferCache() for _ in 1:11])
+end

mwes/mwe_02.jl

@@ -15,12 +15,13 @@ using Interpolations: Extrapolation
 using Parameters
 using Serialization
 using SharedArrays
+using PreallocationTools
 using Pkg
 
 # Export public interface
 export Wing, Section, RamAirWing
 export BodyAerodynamics
-export Solver, solve, solve_base, solve!, VSMSolution
+export Solver, solve, solve_base!, solve!, VSMSolution
 export calculate_results
 export add_section!, set_va!
 export calculate_span, calculate_projected_area

src/VortexStepMethod.jl

@@ -297,16 +297,15 @@ Calculate Aerodynamic Influence Coefficient matrices.
 
 See also: [BodyAerodynamics](@ref), [Model](@ref)
 
-Returns:
-    Tuple of (`AIC_x`, `AIC_y`, `AIC_z`) matrices
+Returns: nothing
 """
-function calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::Model,
+@inline function calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::Model,
                               core_radius_fraction::Float64,
                               va_norm_array::Vector{Float64}, 
                               va_unit_array::Matrix{Float64})
     # Determine evaluation point based on model
-    evaluation_point = model === VSM ? :control_point : :aero_center
-    evaluation_point_on_bound = model === LLT
+    evaluation_point = model == VSM ? :control_point : :aero_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)
@@ -330,13 +329,13 @@ function calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::Model,
                 core_radius_fraction,
                 body_aero.work_vectors
             )
-            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
+                velocity_induced .-= U_2D              
             end
+            body_aero.AIC[:, icp, jring] .= velocity_induced
         end
     end
     return nothing
@@ -347,24 +346,23 @@ end
 
 Calculate circulation distribution for an elliptical wing.
 
-Returns:
-    Vector{Float64}: Circulation distribution along the wing
+Returns: nothing
 """
-function calculate_circulation_distribution_elliptical_wing(body_aero::BodyAerodynamics, gamma_0::Float64=1.0)
+function calculate_circulation_distribution_elliptical_wing(gamma_i, body_aero::BodyAerodynamics, gamma_0::Float64=1.0)
     length(body_aero.wings) == 1 || throw(ArgumentError("Multiple wings not yet implemented"))
 
     wing_span = body_aero.wings[1].span
     @debug "Wing span: $wing_span"
 
     # Calculate y-coordinates of control points
+    # TODO: pre-allocate y
     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)
+    gamma_i .= gamma_0 * sqrt.(1 .- (2 .* y ./ wing_span).^2)
 
     @debug "Calculated circulation distribution: $gamma_i"
-    
-    return gamma_i
+    nothing
 end
 
 """
@@ -413,6 +411,8 @@ function calculate_stall_angle_list(panels::Vector{Panel};
     return stall_angles
 end
 
+const cache_body = [LazyBufferCache() for _ in 1:5]
+
 """
     update_effective_angle_of_attack_if_VSM(body_aero::BodyAerodynamics, gamma::Vector{Float64},
                                           core_radius_fraction::Float64,
@@ -427,7 +427,8 @@ 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(body_aero::BodyAerodynamics, 
+function update_effective_angle_of_attack!(alpha_corrected,
+    body_aero::BodyAerodynamics, 
     gamma::Vector{Float64},
     core_radius_fraction::Float64,
     z_airf_array::Matrix{Float64},
@@ -436,26 +437,53 @@ function update_effective_angle_of_attack_if_VSM(body_aero::BodyAerodynamics,
     va_norm_array::Vector{Float64},
     va_unit_array::Matrix{Float64})
 
-    # Calculate AIC matrices at aerodynamic center using LLT method
-    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 = [
-        AIC_x * gamma,
-        AIC_y * gamma,
-        AIC_z * gamma
-    ]
-    induced_velocity = hcat(induced_velocity...)
+    # Calculate AIC matrices (keep existing optimized view)
+    calculate_AIC_matrices!(body_aero, LLT, core_radius_fraction, va_norm_array, va_unit_array)
+
+    # Get dimensions from existing data
+    n_rows = size(body_aero.AIC, 2)
+    n_cols = size(body_aero.AIC, 3)
+
+    # Preallocate induced velocity array
+    induced_velocity = cache_body[1][va_array]
+
+    # Calculate each component with explicit loops
+    for j in 1:3  # For each x/y/z component
+        for i in 1:n_rows
+            acc = zero(eltype(induced_velocity))  # Type-stable accumulator
+            for k in 1:n_cols
+                acc += body_aero.AIC[j, i, k] * gamma[k]
+            end
+            induced_velocity[i, j] = acc
+        end
+    end
+
+    # In-place relative velocity calculation
+    relative_velocity = cache_body[2][va_array]
+    relative_velocity .= va_array .+ induced_velocity
+
+    # Preallocate and compute dot products manually
+    n = size(relative_velocity, 1)
+    v_normal     = cache_body[3][relative_velocity]
+    v_tangential = cache_body[4][relative_velocity]
 
-    # Calculate relative velocities and angles
-    relative_velocity = va_array + induced_velocity
-    v_normal = sum(z_airf_array .* relative_velocity, dims=2)
-    v_tangential = sum(x_airf_array .* relative_velocity, dims=2)
-    alpha_array = atan.(v_normal ./ v_tangential)
-    return alpha_array
+    @inbounds for i in 1:n
+        vn = 0.0
+        vt = 0.0
+        for j in 1:3
+            vn += z_airf_array[i, j] * relative_velocity[i, j]
+            vt += x_airf_array[i, j] * relative_velocity[i, j]
+        end
+        v_normal[i] = vn
+        v_tangential[i] = vt
+    end
+
+    # Direct angle calculation without temporary arrays
+    @inbounds for i in 1:n
+        alpha_corrected[i] = atan(v_normal[i], v_tangential[i])
+    end
+
+    nothing
 end
 
 """
@@ -501,6 +529,7 @@ function calculate_results(
     cd_array = zeros(n_panels)
     cm_array = zeros(n_panels)
     panel_width_array = zeros(n_panels)
+    alpha_corrected = zeros(n_panels)
 
     # Calculate coefficients for each panel
     for (i, panel) in enumerate(panels)
@@ -515,8 +544,9 @@ function calculate_results(
     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
-        update_effective_angle_of_attack_if_VSM(
+    if aerodynamic_model_type === VSM
+        update_effective_angle_of_attack!(
+            alpha_corrected,
             body_aero,
             gamma_new,
             core_radius_fraction,
@@ -526,10 +556,8 @@ function calculate_results(
             va_norm_array,
             va_unit_array
         )
-    elseif aerodynamic_model_type === LLT
-        alpha_array
-    else
-        throw(ArgumentError("Unknown aerodynamic model type, should be LLT or VSM"))
+    elseif aerodynamic_model_type == LLT
+        alpha_corrected .= alpha_array
     end
 
     # Verify va is not distributed