Merge branch 'main' into CI

Author: ufechner7

Project.toml

@@ -4,6 +4,7 @@ authors = ["1-Bart-1 <bart@vandelint.net>"]
 version = "0.1.0"
 
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"

src/filament.jl

@@ -35,35 +35,45 @@ function velocity_3D_bound_vortex!(
     filament::BoundFilament,
     XVP::Vector{Float64},
     gamma::Float64,
-    core_radius_fraction::Float64
+    core_radius_fraction::Float64,
+    work_vectors::NTuple{10, Vector{Float64}}
 )
+    r1, r2, r1Xr2, r1Xr0, r2Xr0, r1r2norm, r1_proj, r2_proj, r1_projXr2_proj, vel_ind_proj = work_vectors
     r0 = filament.r0
-    r1 = XVP - filament.x1
-    r2 = XVP - filament.x2
+    r1 .= XVP .- filament.x1
+    r2 .= XVP .- filament.x2
 
     # Cut-off radius
     epsilon = core_radius_fraction * norm(r0)
+
+    cross3!(r1Xr2, r1, r2)
+    cross3!(r1Xr0, r1, r0)
+    r1r2norm .= r1./norm(r1) .- r2./norm(r2)
 
     # Check point location relative to filament
-    if norm(cross(r1, r0)) / norm(r0) > epsilon
-        vel .= (gamma / (4π)) * cross(r1, r2) / (norm(cross(r1, r2))^2) * 
-               dot(r0, r1/norm(r1) - r2/norm(r2))
-    elseif norm(cross(r1, r0)) / norm(r0) == 0
+    if norm(r1Xr0) / norm(r0) > epsilon
+        vel .= (gamma / (4π)) .* r1Xr2 ./ (norm(r1Xr2)^2) .* 
+            dot(r0, r1r2norm)
+    elseif norm(r1Xr0) / norm(r0) == 0
         vel .= zeros(3)
     else
         @info "inside core radius"
-        @info "distance from control point to filament: $(norm(cross(r1, r0)) / norm(r0))"
+        @info "distance from control point to filament: $(norm(r1Xr0) / norm(r0))"
 
         # Project onto core radius
-        r1_proj = dot(r1, r0) * r0 / (norm(r0)^2) + 
-                  epsilon * cross(r1, r0) / norm(cross(r1, r0))
-        r2_proj = dot(r2, r0) * r0 / (norm(r0)^2) + 
-                  epsilon * cross(r2, r0) / norm(cross(r2, r0))
-        
-        vel_ind_proj = (gamma / (4π)) * cross(r1_proj, r2_proj) / (norm(cross(r1_proj, r2_proj))^2) * 
-                      dot(r0, r1_proj/norm(r1_proj) - r2_proj/norm(r2_proj))
-        
-        vel .= norm(cross(r1, r0)) / (norm(r0) * epsilon) * vel_ind_proj
+        cross3!(r2Xr0, r2, r0)
+        r1_proj .= dot(r1, r0) .* r0 ./ (norm(r0)^2) .+ 
+                  epsilon .* r1Xr0 ./ norm(r1Xr0)
+        r2_proj .= dot(r2, r0) .* r0 ./ (norm(r0)^2) .+ 
+                  epsilon .* r2Xr0 ./ norm(r2Xr0)
+        cross3!(r1_projXr2_proj, r1_proj, r2_proj)
+
+        vel_ind_proj .= (gamma / (4π)) .* r1_projXr2_proj ./ (norm(r1_projXr2_proj)^2) .* 
+                      dot(r0, r1_proj/norm(r1_proj) .- r2_proj/norm(r2_proj))
+
+        @show vel_ind_proj r1_projXr2_proj
+
+        vel .= norm(r1Xr0) ./ (norm(r0) * epsilon) .* vel_ind_proj
     end
     nothing
 end

src/panel.jl

@@ -374,7 +374,8 @@ function calculate_velocity_induced_single_ring_semiinfinite(
                     filament,
                     evaluation_point,
                     gamma,
-                    core_radius_fraction
+                    core_radius_fraction,
+                    work_vectors
                 )
             end
         elseif i ∈ (2, 3)  # trailing filaments

test/runtests.jl

@@ -4,4 +4,8 @@ cd("..")
 println("Running tests...")
 @testset verbose = true "Testing VortexStepMethod..." begin
     include("test_bound_filament.jl")
+    include("test_panel.jl")
+    include("test_semi_infinite_filament.jl")
+    include("test_wing_aerodynamics.jl")
+    include("test_wing_geometry.jl")
 end

test/test_bound_filament.jl

@@ -1,5 +1,6 @@
 using VortexStepMethod: BoundFilament, velocity_3D_bound_vortex!
 using LinearAlgebra
+using BenchmarkTools: @benchmarkable
 using Test
 
 # Test helper functions
@@ -27,6 +28,25 @@ end
     gamma = 1.0
     core_radius_fraction = 0.01
     induced_velocity = zeros(3)
+    work_vectors = ntuple(_ -> Vector{Float64}(undef, 3), 10)
+
+    @testset "Non-allocating" begin
+        filament = create_test_filament()
+        control_point = [0.5, 1.0, 0.0]
+        induced_velocity = zeros(3)
+        
+        b = @benchmarkable velocity_3D_bound_vortex!(
+            $induced_velocity,
+            $filament,
+            $control_point,
+            $gamma,
+            $core_radius_fraction,
+            $work_vectors
+        )
+        result = run(b)
+        @test result.allocs == 0
+        @test result.memory == 0
+    end
 
     @testset "Calculate Induced Velocity" begin
         filament = create_test_filament()
@@ -37,7 +57,8 @@ end
             filament,
             control_point,
             gamma,
-            core_radius_fraction
+            core_radius_fraction,
+            work_vectors
         )
         analytical_velocity = analytical_solution(control_point, gamma)
 
@@ -58,7 +79,8 @@ end
                 filament,
                 point,
                 gamma,
-                core_radius_fraction
+                core_radius_fraction,
+                work_vectors
             )
             @test all(isapprox.(induced_velocity, zeros(3), atol=1e-10))
             @test !any(isnan.(induced_velocity))
@@ -74,9 +96,9 @@ end
             filament,
             control_point,
             gamma,
-            core_radius_fraction
+            core_radius_fraction,
+            work_vectors
         )
-        
         @test !any(isnan.(induced_velocity))
         @test isapprox(induced_velocity[1], 0.0, atol=1e-8)
         @test abs(induced_velocity[2]) < 1e-8
@@ -92,7 +114,8 @@ end
             filament,
             control_point,
             gamma,
-            core_radius_fraction
+            core_radius_fraction,
+            work_vectors
         )
 
         @test !any(isnan.(induced_velocity))
@@ -105,9 +128,9 @@ end
         filament = create_test_filament()
         control_point = [0.5, 1.0, 0.0]
 
-        velocity_3D_bound_vortex!(v1, filament, control_point, 1.0, core_radius_fraction)
-        velocity_3D_bound_vortex!(v2, filament, control_point, 2.0, core_radius_fraction)
-        velocity_3D_bound_vortex!(v4, filament, control_point, 4.0, core_radius_fraction)
+        velocity_3D_bound_vortex!(v1, filament, control_point, 1.0, core_radius_fraction, work_vectors)
+        velocity_3D_bound_vortex!(v2, filament, control_point, 2.0, core_radius_fraction, work_vectors)
+        velocity_3D_bound_vortex!(v4, filament, control_point, 4.0, core_radius_fraction, work_vectors)
 
         @test isapprox(v4, 2 * v2)
         @test isapprox(v4, 4 * v1)
@@ -116,8 +139,8 @@ end
     @testset "Symmetry" begin
         filament = BoundFilament([-1.0, 0.0, 0.0], [1.0, 0.0, 0.0])
 
-        velocity_3D_bound_vortex!(v1, filament, [0.0, 1.0, 0.0], gamma, core_radius_fraction)
-        velocity_3D_bound_vortex!(v2, filament, [0.0, -1.0, 0.0], gamma, core_radius_fraction)
+        velocity_3D_bound_vortex!(v1, filament, [0.0, 1.0, 0.0], gamma, core_radius_fraction, work_vectors)
+        velocity_3D_bound_vortex!(v2, filament, [0.0, -1.0, 0.0], gamma, core_radius_fraction, work_vectors)
 
         @test isapprox(v1, -v2)
     end
@@ -134,7 +157,7 @@ end
 
         velocities = [zeros(3) for p in points]
         [
-            velocity_3D_bound_vortex!(velocities[i], filament, p, gamma, core_radius_fraction)
+            velocity_3D_bound_vortex!(velocities[i], filament, p, gamma, core_radius_fraction, work_vectors)
             for (i, p) in enumerate(points)
         ]
 
@@ -147,7 +170,7 @@ end
         @test norm(velocities[2]) > norm(velocities[3])
 
         # Check continuity around core radius
-        @test_broken isapprox(velocities[1], velocities[2], rtol=1e-3)
+        @test isapprox(velocities[1], velocities[2], rtol=1e-2)
 
         # Check non-zero velocities
         @test !all(isapprox.(velocities[1], zeros(3), atol=1e-10))
@@ -161,7 +184,8 @@ end
             filament,
             [0.5, -core_radius_fraction, 0.0],
             gamma,
-            core_radius_fraction
+            core_radius_fraction,
+            work_vectors
         )
         @test isapprox(velocities[2], -v_neg)
     end

test/test_wing_aerodynamics.jl

@@ -63,7 +63,7 @@ include("utils.jl")
         cd_cm = calculate_cd_cm(panel, alpha[i])
         drag[i] = dyn_visc * cd_cm[1] * panel.chord
         moment[i] = dyn_visc * cd_cm[2] * panel.chord^2
-        @info "lift: $lift, drag: $drag, moment: $moment"
+        # @info "lift: $lift, drag: $drag, moment: $moment"
     end
     Fmag = hcat(lift, drag, moment)
 
@@ -104,7 +104,6 @@ include("utils.jl")
     @test length(results_NEW["cd_distribution"]) == length(wing_aero.panels)
 end
 
-
 @testset "Induction Matrix Creation" begin
     # Setup
     n_panels = 3
@@ -206,6 +205,149 @@ end
         @test isapprox(MatrixU, AIC_x, atol=1e-8)
         @test isapprox(MatrixV, -AIC_y, atol=1e-8)
         @test isapprox(MatrixW, AIC_z, atol=1e-8)
-        @show MatrixU MatrixV MatrixW
     end
-end
+end
+
+
+@testset "Wing Geometry Creation" begin
+    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)
+        @debug "AR: $AR"
+        Umag = 20.0
+        aoa = 5.7106 * π / 180
+        Uinf = [cos(aoa), 0.0, sin(aoa)] .* Umag
+    
+        coord = if wing_type == "rectangular"
+            twist = range(-0.5, 0.5, length=N)
+            beta = range(-2, 2, length=N)
+            generate_coordinates_rect_wing(
+                fill(max_chord, N),
+                span,
+                twist,
+                beta,
+                N,
+                "lin"
+            )
+        elseif wing_type == "curved"
+            generate_coordinates_curved_wing(
+                max_chord, span, π/4, 5, N, "cos"
+            )
+        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")
+        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"
+            )
+        end
+        wing_aero = WingAerodynamics([wing])
+        set_va!(wing_aero, (Uinf, 0.0))
+        
+        return wing_aero, coord, Uinf, model
+    end
+
+    for model in ["VSM", "LLT"]
+        @debug "model: $model"
+        for wing_type in ["rectangular", "curved", "elliptical"]
+            @debug "wing_type: $wing_type"
+            wing_aero, coord, Uinf, model = create_geometry(
+                model=model, wing_type=wing_type
+            )
+            
+            # Generate geometry
+            expected_controlpoints, expected_rings, expected_bladepanels, 
+                expected_ringvec, expected_coord_L = create_geometry_general(
+                coord, Uinf, div(size(coord,1), 2), "5fil", model
+            )
+
+            for i in 1:length(wing_aero.panels)
+                @debug "i: $i"
+                # Handle control points
+                index_reversed = length(wing_aero.panels) - i + 1
+                panel = wing_aero.panels[index_reversed]
+                
+                evaluation_point = if model == "VSM"
+                    panel.control_point
+                else  # LLT
+                    panel.aerodynamic_center
+                end
+
+                @test isapprox(evaluation_point, expected_controlpoints[i]["coordinates"], atol=1e-4)
+                @test isapprox(panel.chord, expected_controlpoints[i]["chord"], atol=1e-4)
+                @test isapprox(panel.x_airf, expected_controlpoints[i]["normal"], atol=1e-4)
+                @test isapprox(panel.y_airf, expected_controlpoints[i]["tangential"], atol=1e-4)
+                @test isapprox(
+                    hcat(panel.x_airf, panel.y_airf, panel.z_airf),
+                    expected_controlpoints[i]["airf_coord"],
+                    atol=1e-4
+                )
+                
+                if model == "VSM"
+                    @test isapprox(
+                        panel.aerodynamic_center,
+                        expected_controlpoints[i]["coordinates_aoa"],
+                        atol=1e-4
+                    )
+                end
+
+                # Handle rings
+                expected_ring_i = expected_rings[i]
+                expected_ring_i_list = [
+                    expected_ring_i[1],
+                    expected_ring_i[2],
+                    expected_ring_i[3],
+                    expected_ring_i[4],
+                    expected_ring_i[5]
+                ]
+
+                filaments = panel.filaments
+                filament_list = [
+                    filaments[1],
+                    filaments[3],
+                    filaments[5],
+                    filaments[2],
+                    filaments[4]
+                ]
+
+                for (j, fil) in enumerate(filament_list)
+                    if j == 1  # bound filaments
+                        @test isapprox(fil.x1, expected_ring_i_list[j]["x1"], atol=1e-4)
+                        @test isapprox(fil.x2, expected_ring_i_list[j]["x2"], atol=1e-4)
+                    elseif j ∈ (2, 4)  # trailing filaments
+                        @test isapprox(fil.x1, expected_ring_i_list[j]["x1"], atol=1e-4)
+                        @test isapprox(fil.x2, expected_ring_i_list[j]["x2"], atol=1e-4)
+                    else  # semi-infinite filaments
+                        @test isapprox(fil.x1, expected_ring_i_list[j]["x1"], atol=1e-4)
+                    end
+                end
+
+                # Handle bladepanels
+                exp_bladepanels = expected_bladepanels[i]
+                @test isapprox(panel.LE_point_2, exp_bladepanels["p1"], atol=1e-4)
+                @test isapprox(panel.LE_point_1, exp_bladepanels["p2"], atol=1e-4)
+                @test isapprox(panel.TE_point_1, exp_bladepanels["p3"], atol=1e-4)
+                @test isapprox(panel.TE_point_2, exp_bladepanels["p4"], atol=1e-4)
+
+                # Handle ringvec
+                exp_ringvec = expected_ringvec[i]
+                r0 = panel.bound_point_1 - panel.bound_point_2
+                r3 = evaluation_point - (panel.bound_point_1 + panel.bound_point_2) / 2
+                @test isapprox(r0, exp_ringvec["r0"], atol=1e-4)
+                @test isapprox(r3, exp_ringvec["r3"], atol=1e-4)
+
+                # Handle coord_L
+                @test all(isapprox.(panel.aerodynamic_center, expected_coord_L[:, i]))
+            end
+        end
+    end
+end