Merge branch 'main' into lin

Author: 1-Bart-1

docs/make.jl

@@ -28,6 +28,7 @@ makedocs(;
         "Private Functions" => "private_functions.md",
         "Private Types" => "private_types.md",
         "Reference Frames" => "reference_frames.md",
+        "Tips and tricks" => "tips_and_tricks.md",
         "Glossary" => "glossary.md"
     ],
 )

docs/src/tips_and_tricks.md

@@ -0,0 +1,29 @@
+## Tips and tricks
+
+### What can this model simulate
+The following bodies can be simulated:
+
+- conventional bodies, consisting of one or more wings
+- leading edge inflatable (LEI) kites
+- RAM-air kites
+
+To build the geometry of a RAM-air kite, a 3D .obj file can be used as input. In addition a `.dat` file is needed.
+It should have two columns, one for the `x` and one for the `y` coordinate of the 2D polar that is used.
+
+### RAM-air kite model
+If running the example `ram_air_kite.jl` fails, try to run the `cleanup.jl` script and then try again. Background: this example caches the calculated polars. Reading cached polars can fail after an update.
+
+### Output formats
+Currently, the `solve!()` function returns the results as [VSMSolution](@ref) struct. The function solve() returns a
+dictionary with the results. The `solve!()` function is faster, and the `solve()` contains many more entries, therefore
+the first function is good for integration in dynamic models and the second one better suited for aerodynamic analysis.
+
+### Performance
+Calling `init!(body_aero; init_aero=false)` is very fast. After calling `deform!(wing)`, you have to run `init!(body_aero; init_aero=false)` to apply the deformed wing to the body aerodynamics. This is in turn necessary for the linearization from deformation to aerodynamic coefficients for RAM-air kites.
+
+### Building the documentation locally
+You can build the documentation locally after checking out the source code with git, launching Julia and executing:
+```
+include("scripts/build_docu.jl")
+```
+A browser window should pop up automatically.

src/solver.jl

@@ -5,30 +5,39 @@
 Struct for storing the solution of the [solve!](@ref) function. Must contain all info needed by `KiteModels.jl`.
 
 # Attributes
-- gamma_distribution::Union{Nothing, Vector{Float64}}: Vector containing the panel circulations
+- `panel_width_array`::Vector{Float64}: Width of the panels [m]
+- cl_array::Vector{Float64}: Lift coefficients of the panels [-]
+- cd_array::Vector{Float64}: Drag coefficients of the panels [-]
+- cm_array::Vector{Float64}: Pitching moment coefficients of the panels [-]
+- panel_lift::Vector{Float64}: Lift force of the panels [N]
+- panel_drag::Vector{Float64}: Drag force of the panels [N]
+- panel_moment::Vector{Float64}: Pitching moment around the spanwise vector of the panels [Nm]
+- `f_body_3D`::Matrix{Float64}: Matrix of the aerodynamic forces (x, y, z vectors) [N]
+- `m_body_3D`::Matrix{Float64}: Matrix of the aerodynamic moments [Nm]
+- `gamma_distribution`::Union{Nothing, Vector{Float64}}: Vector containing the panel circulations.
 - force::MVec3: Aerodynamic force vector in KB reference frame [N]
 - moment::MVec3: Aerodynamic moments [Mx, My, Mz] around the reference point [Nm]
 - force_coefficients::MVec3: Aerodynamic force coefficients [CFx, CFy, CFz] [-]
-- moment_coefficients::MVec3: Aerodynamic moment coefficients [CMx, CMy, CMz] [-]
-- moment_distribution::Vector{Float64}: Pitching moments around the spanwise vector of each panel. [Nm]
-- moment_coeff_dist::Vector{Float64}: Pitching moment coefficient around the spanwise vector of each panel. [-]
-- solver_status::SolverStatus: enum, see [SolverStatus](@ref)
+- `moment_coefficients`::MVec3: Aerodynamic moment coefficients [CMx, CMy, CMz] [-]
+- `moment_distribution`::Vector{Float64}: Pitching moments around the spanwise vector of each panel. [Nm]
+- `moment_coefficient_distribution`::Vector{Float64}: Pitching moment coefficient around the spanwise vector of each panel. [-]
+- `solver_status`::SolverStatus: enum, see [SolverStatus](@ref)
 """
 @with_kw mutable struct VSMSolution{P}
     ### private vectors of solve_base!
-    x_airf_array::Matrix{Float64} = zeros(P, 3)
-    y_airf_array::Matrix{Float64} = zeros(P, 3)
-    z_airf_array::Matrix{Float64} = zeros(P, 3)
-    va_array::Matrix{Float64} = zeros(P, 3)
-    chord_array::MVector{P, Float64} = zeros(MVector{P, Float64})
-    ###
-    panel_width_array::MVector{P, Float64} = zeros(MVector{P, Float64})
-    cl_array::MVector{P, Float64} = zeros(MVector{P, Float64})
-    cd_array::MVector{P, Float64} = zeros(MVector{P, Float64})
-    cm_array::MVector{P, Float64} = zeros(MVector{P, Float64})
-    panel_lift::Matrix{Float64} = zeros(P,1)
-    panel_drag::Matrix{Float64} = zeros(P,1)
-    panel_moment::Matrix{Float64} = zeros(P,1)
+    _x_airf_array::Matrix{Float64} = zeros(P, 3)
+    _y_airf_array::Matrix{Float64} = zeros(P, 3)
+    _z_airf_array::Matrix{Float64} = zeros(P, 3)
+    _va_array::Matrix{Float64} = zeros(P, 3)
+    _chord_array::Vector{Float64} = zeros(P)
+    ### end of private vectors
+    panel_width_array::Vector{Float64} = zeros(P)
+    cl_array::Vector{Float64} = zeros(P)
+    cd_array::Vector{Float64} = zeros(P)
+    cm_array::Vector{Float64} = zeros(P)
+    panel_lift::Vector{Float64} = zeros(P)
+    panel_drag::Vector{Float64} = zeros(P)
+    panel_moment::Vector{Float64} = zeros(P)
     f_body_3D::Matrix{Float64} = zeros(3, P)
     m_body_3D::Matrix{Float64} = zeros(3, P)
     gamma_distribution::Union{Nothing, Vector{Float64}} = nothing
@@ -196,9 +205,9 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=
     panel_moment = solver.sol.panel_moment
 
     # Compute using fused broadcasting (no intermediate allocations)
-    @. lift = cl_array * 0.5 * density * v_a_array^2 * solver.sol.chord_array
-    @. drag = cd_array * 0.5 * density * v_a_array^2 * solver.sol.chord_array
-    @. panel_moment = cm_array * 0.5 * density * v_a_array^2 * solver.sol.chord_array
+    @. lift = cl_array * 0.5 * density * v_a_array^2 * solver.sol._chord_array
+    @. drag = cd_array * 0.5 * density * v_a_array^2 * solver.sol._chord_array
+    @. panel_moment = cm_array * 0.5 * density * v_a_array^2 * solver.sol._chord_array
 
     # Calculate alpha corrections based on model type
     if aerodynamic_model_type == VSM                             # 64 bytes
@@ -207,9 +216,9 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=
             body_aero,
             gamma_new,
             solver.core_radius_fraction,
-            solver.sol.z_airf_array,
-            solver.sol.x_airf_array,
-            solver.sol.va_array,
+            solver.sol._z_airf_array,
+            solver.sol._x_airf_array,
+            solver.sol._va_array,
             solver.br.va_norm_array,
             solver.br.va_unit_array
         )
@@ -334,11 +343,11 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n
         solver.mu,
         solver.lr.alpha_array,
         solver.lr.v_a_array,
-        solver.sol.chord_array,
-        solver.sol.x_airf_array,
-        solver.sol.y_airf_array,
-        solver.sol.z_airf_array,
-        solver.sol.va_array,
+        solver.sol._chord_array,
+        solver.sol._x_airf_array,
+        solver.sol._y_airf_array,
+        solver.sol._z_airf_array,
+        solver.sol._va_array,
         solver.br.va_norm_array,
         solver.br.va_unit_array,
         body_aero.panels,
@@ -365,31 +374,31 @@ function solve_base!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribu
     relaxation_factor = solver.relaxation_factor
 
     # Clear arrays
-    solver.sol.x_airf_array .= 0
-    solver.sol.y_airf_array .= 0
-    solver.sol.z_airf_array .= 0
-    solver.sol.va_array .= 0
-    solver.sol.chord_array .= 0
+    solver.sol._x_airf_array .= 0
+    solver.sol._y_airf_array .= 0
+    solver.sol._z_airf_array .= 0
+    solver.sol._va_array .= 0
+    solver.sol._chord_array .= 0
 
     # Fill arrays from panels
     for (i, panel) in enumerate(panels)
-        solver.sol.x_airf_array[i, :] .= panel.x_airf
-        solver.sol.y_airf_array[i, :] .= panel.y_airf
-        solver.sol.z_airf_array[i, :] .= panel.z_airf
-        solver.sol.va_array[i, :] .= panel.va
-        solver.sol.chord_array[i] = panel.chord
+        solver.sol._x_airf_array[i, :] .= panel.x_airf
+        solver.sol._y_airf_array[i, :] .= panel.y_airf
+        solver.sol._z_airf_array[i, :] .= panel.z_airf
+        solver.sol._va_array[i, :] .= panel.va
+        solver.sol._chord_array[i] = panel.chord
     end
 
     # Calculate unit vectors
-    calc_norm_array!(solver.br.va_norm_array, solver.sol.va_array)
-    solver.br.va_unit_array .= solver.sol.va_array ./ solver.br.va_norm_array
+    calc_norm_array!(solver.br.va_norm_array, solver.sol._va_array)
+    solver.br.va_unit_array .= solver.sol._va_array ./ solver.br.va_norm_array
 
     # Calculate AIC matrices
     calculate_AIC_matrices!(body_aero, solver.aerodynamic_model_type, solver.core_radius_fraction, solver.br.va_norm_array,
                             solver.br.va_unit_array)
 
     # Initialize gamma distribution
-    gamma_initial = solver.cache_base[1][solver.sol.chord_array]
+    gamma_initial = solver.cache_base[1][solver.sol._chord_array]
     if isnothing(gamma_distribution)
         if solver.type_initial_gamma_distribution == ELLIPTIC
             calculate_circulation_distribution_elliptical_wing(gamma_initial, body_aero)
@@ -430,11 +439,11 @@ function gamma_loop!(
     relaxation_factor::Float64;
     log::Bool = true
 )
-    va_array = solver.sol.va_array
-    chord_array = solver.sol.chord_array
-    x_airf_array = solver.sol.x_airf_array
-    y_airf_array = solver.sol.y_airf_array
-    z_airf_array = solver.sol.z_airf_array
+    va_array = solver.sol._va_array
+    chord_array = solver.sol._chord_array
+    x_airf_array = solver.sol._x_airf_array
+    y_airf_array = solver.sol._y_airf_array
+    z_airf_array = solver.sol._z_airf_array
     solver.lr.converged   = false
     n_panels    = length(body_aero.panels)
     solver.lr.alpha_array .= body_aero.alpha_array

test/bench.jl

@@ -129,11 +129,11 @@ using LinearAlgebra
                 solver = Solver{P}(
                     aerodynamic_model_type=model
                 )
-                solver.sol.va_array .= va_array
-                solver.sol.chord_array .= chord_array
-                solver.sol.x_airf_array .= x_airf_array
-                solver.sol.y_airf_array .= y_airf_array
-                solver.sol.z_airf_array .= z_airf_array
+                solver.sol._va_array .= va_array
+                solver.sol._chord_array .= chord_array
+                solver.sol._x_airf_array .= x_airf_array
+                solver.sol._y_airf_array .= y_airf_array
+                solver.sol._z_airf_array .= z_airf_array
                 result = @benchmark gamma_loop!(
                     $solver,
                     $body_aero,

test/bench2.jl

@@ -0,0 +1,56 @@
+using LinearAlgebra
+using VortexStepMethod
+using BenchmarkTools
+using Test
+
+using Pkg
+
+if !("CSV" ∈ keys(Pkg.project().dependencies))
+    using TestEnv
+    TestEnv.activate()
+end
+
+# 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_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)
+
+# Step 4: Initialize solvers for both LLT and VSM methods
+P = length(wa.panels)
+vsm_solver = Solver{P}(aerodynamic_model_type=VSM)
+
+@testset "Allocation Tests for solve() and solve!()" begin
+    # Step 5: Solve using both methods
+    result = @benchmark  solve_base!($vsm_solver, $wa, nothing)  # 51 allocations
+    @test result.allocs <= 55
+    # time Python: 32.0 ms  Ryzen 7950x
+    # time Julia:   0.45 ms Ryzen 7950x
+    result = @benchmark  sol = solve!($vsm_solver, $wa, nothing) # 85 allocations
+    @test result.allocs <= 89
+end
+nothing
+sol = solve!(vsm_solver, wa, nothing) # 85 allocations