From 9879fc08ec3289e0e637942b2b6888c176e53f69 Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 17 Aug 2026 00:07:10 +0200 Subject: [PATCH 1/7] Store AIC component-last so the induced-velocity products reach BLAS MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `BodyAerodynamics.AIC` was `(3, n_panels, n_panels)`, so each component slice `AIC[k, :, :]` had `stride1 == 3`. That is not a BLAS layout, so the three `mul!` calls in the gamma loop silently fell back to the generic Julia matmul — over 247 iterations that was ~80% of `solve!`. Storing it as `(n_panels, n_panels, 3)` makes every slice contiguous and the products dispatch to `gemv`: before after speedup INVISCID n=60 5.92 ms 1.72 ms 3.43x INVISCID n=120 26.10 ms 4.89 ms 5.34x POLAR_VECTORS n=60 6.44 ms 2.11 ms 3.06x POLAR_VECTORS n=120 21.06 ms 6.05 ms 3.48x Also swaps the AIC build loop nest to `jring` outer so the writes are unit-stride and the filament tuple hoists, replaces the hand-rolled triple loop in `update_effective_angle_of_attack!` with three `mul!` calls, and drops the AIC-as-shape-template for `f_body_3D`/`m_body_3D`, which no longer has a 3xn slice to borrow. Allocation counts are unchanged and `calc_forces! is zero-alloc` still holds. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 5 +++ src/body_aerodynamics.jl | 45 +++++++------------ src/solver.jl | 6 +-- .../test_body_aerodynamics.jl | 4 +- 4 files changed, 27 insertions(+), 33 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f2737cb6..ba4f401c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,11 @@ changes format without re-running the airfoil solver that produced it. ### Changed +- BREAKING: `BodyAerodynamics.AIC` is stored as `(n_panels, n_panels, 3)` instead of + `(3, n_panels, n_panels)`, so each component slice `AIC[:, :, k]` is contiguous and + the induced-velocity products reach BLAS `gemv` instead of the generic fallback. + `solve!` is 3.1–5.3× faster (n=120, VSM: 26.1 ms → 4.9 ms inviscid, 21.1 ms → + 6.0 ms with polars). Code reading `AIC[k, i, j]` must become `AIC[i, j, k]`. - `read_node_table` parses into a preallocated matrix instead of `reduce(vcat, …)` over a generator, which was quadratic in the row count: ~21× faster on a 16 MB surface table (2.49 s → 0.12 s), benefiting every existing dataset. diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 044d4df7..87b514ec 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -16,7 +16,8 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru - `alpha_dist::MVector{P, Float64}` = zeros(Float64, P) - `v_a_dist::MVector{P, Float64}` = zeros(Float64, P) - `work_vectors`::NTuple{10, MVec3} = ntuple(_ -> zeros(MVec3), 10) -- `AIC::Array{Float64, 3}` = zeros(3, P, P) +- `AIC::Array{Float64, 3}` = zeros(P, P, 3): influence coefficients, component last so + that each `AIC[:, :, k]` slice is a contiguous BLAS matrix - `projected_area::Float64` = 1.0: The area projected onto the xy-plane of the kite body reference frame [m²] - `c_ref::Float64` = 1.0: Reference chord length (max panel chord) [m] - `y::MVector{P, Float64}` = MVector{P,Float64}(zeros(P)) @@ -35,7 +36,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru alpha_dist::MVector{P, T} = zeros(MVector{P, T}) v_a_dist::MVector{P, T} = zeros(MVector{P, T}) work_vectors::NTuple{10, MVector{3, T}} = ntuple(_ -> zeros(MVector{3, T}), 10) - AIC::Array{T, 3} = zeros(T, 3, P, P) + AIC::Array{T, 3} = zeros(T, P, P, 3) projected_area::T = one(T) c_ref::T = one(T) y::MVector{P, T} = zeros(MVector{P, T}) @@ -397,12 +398,13 @@ Returns: nothing va_norm = wake_speed # Calculate influence coefficients - for icp in eachindex(body_aero.panels) - panel_icp = body_aero.panels[icp] - ep = evaluation_point == :control_point ? panel_icp.control_point : panel_icp.aero_center - for jring in eachindex(body_aero.panels) - panel_jring = body_aero.panels[jring] - filaments = panel_jring.filaments + for jring in eachindex(body_aero.panels) + panel_jring = body_aero.panels[jring] + filaments = panel_jring.filaments + for icp in eachindex(body_aero.panels) + panel_icp = body_aero.panels[icp] + ep = evaluation_point == :control_point ? panel_icp.control_point : + panel_icp.aero_center calculate_velocity_induced_single_ring_semiinfinite!( velocity_induced, tempvel, @@ -421,7 +423,9 @@ Returns: nothing calculate_velocity_induced_bound_2D!(U_2D, panel_jring, ep, body_aero.work_vectors) velocity_induced .-= U_2D end - body_aero.AIC[:, icp, jring] .= velocity_induced + @inbounds for k in 1:3 + body_aero.AIC[icp, jring, k] = velocity_induced[k] + end end end return nothing @@ -478,25 +482,11 @@ function update_effective_angle_of_attack!(alpha_corrected, va_norm_array, va_unit_array) - # 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 = body_aero.cache[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 + for k in 1:3 + mul!(view(induced_velocity, :, k), view(body_aero.AIC, :, :, k), gamma) end # In-place relative velocity calculation @@ -762,9 +752,8 @@ function calculate_results( cl_prescribed_va = body_aero.cache[10][alpha_dist] cd_prescribed_va = body_aero.cache[11][alpha_dist] cs_prescribed_va = body_aero.cache[12][alpha_dist] - panel_view_3xn = @view body_aero.AIC[:, :, 1] - f_body_3D = body_aero.cache[13][panel_view_3xn] - m_body_3D = body_aero.cache[14][panel_view_3xn] + f_body_3D = body_aero.cache[13][alpha_dist, (3, length(alpha_dist))] + m_body_3D = body_aero.cache[14][alpha_dist, (3, length(alpha_dist))] alpha_geometric = body_aero.cache[15][alpha_dist] fill!(f_body_3D, 0.0) diff --git a/src/solver.jl b/src/solver.jl index 9d5ff060..d008e544 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -904,9 +904,9 @@ function gamma_loop!( v_normal_array = solver.cache[10][solver.lr.gamma_new] v_tangential_array = solver.cache[11][solver.lr.gamma_new] - AIC_x = @view body_aero.AIC[1, :, :] - AIC_y = @view body_aero.AIC[2, :, :] - AIC_z = @view body_aero.AIC[3, :, :] + AIC_x = @view body_aero.AIC[:, :, 1] + AIC_y = @view body_aero.AIC[:, :, 2] + AIC_z = @view body_aero.AIC[:, :, 3] velocity_view_x = @view induced_velocity_all[:, 1] velocity_view_y = @view induced_velocity_all[:, 2] diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl index 7ecfd6f8..ab705b2b 100644 --- a/test/body_aerodynamics/test_body_aerodynamics.jl +++ b/test/body_aerodynamics/test_body_aerodynamics.jl @@ -72,7 +72,7 @@ end 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, :, :] + 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) @@ -107,7 +107,7 @@ end va_norm_array, va_unit_array ) - AIC_x, AIC_y, AIC_z = body_aero.AIC[1, :, :], body_aero.AIC[2, :, :], body_aero.AIC[3, :, :] + 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) From ae73645f67b85d471e2f711e296dd9fb8b6b87a7 Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 17 Aug 2026 00:38:15 +0200 Subject: [PATCH 2/7] Drop redundant work from the filament velocity kernels MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The three induced-velocity kernels run n^2 times per `solve!`, and each was recomputing quantities it either already had or did not need: - `r0` and `|r0|` are `reinit!`-maintained on every `BoundFilament`, but `velocity_3D_trailing_vortex!` rebuilt both from `x1`/`x2` each call and `velocity_3D_bound_vortex!` re-normed the stored vector. - The core-radius cutoff only needs `|r_perp|`, and `r_perp` is a projection onto `r0` (or onto `Vf`), so its length is `|r1.r0|/|r0|` — no need to form the vector. It is still built inside the core branch, where it is used. - `r1 x r2`, the `r1/|r1| - r2/|r2|` difference and `r2 x r0` were computed before the branch that selects between the regular and core-radius forms, but each is read on only one side of it. Bit-identical output: max deviation 0.0 over 28800 calls against the previous implementation, and the wing coefficients are unchanged. INVISCID n=60 1724.0 us -> 1173.3 us 1.47x INVISCID n=120 4889.0 us -> 3277.7 us 1.49x POLAR_VECTORS n=60 2106.2 us -> 1428.4 us 1.47x POLAR_VECTORS n=120 6047.6 us -> 3908.3 us 1.55x The kernels keep their `MVector` scratch. An `SVector` rewrite is 1.14x faster again for `Float64` but 1.2-1.4x *slower* for `ForwardDiff.Dual`, which `linearize` depends on, so it is not worth the trade. Also fixes two stale defaults in the `Solver` docstring: `type_initial_gamma_distribution` is `ZEROS`, not `ELLIPTIC`, and `core_radius_fraction` is `0.05`, not `1e-20`. Co-Authored-By: Claude Opus 5 (1M context) --- src/filament.jl | 55 +++++++++++++++++++++---------------------------- src/solver.jl | 4 ++-- 2 files changed, 25 insertions(+), 34 deletions(-) diff --git a/src/filament.jl b/src/filament.jl index c1a5d5d5..6ea8bceb 100644 --- a/src/filament.jl +++ b/src/filament.jl @@ -55,24 +55,23 @@ function velocity_3D_bound_vortex!( r1, r2, r1Xr2, r1Xr0, r2Xr0, r1r2norm, r1_proj, r2_proj, r1_projXr2_proj, vel_ind_proj = work_vectors r0 = filament.r0 + nr0 = filament.length r1 .= XVP .- filament.x1 r2 .= XVP .- filament.x2 - # Cut-off radius - nr0 = norm3(r0) epsilon = core_radius_fraction * nr0 - cross3!(r1Xr2, r1, r2) cross3!(r1Xr0, r1, r0) - nr1 = norm3(r1) - nr2 = norm3(r2) - @inbounds for k in 1:3 - r1r2norm[k] = r1[k]/nr1 - r2[k]/nr2 - end # Check point location relative to filament nr1Xr0 = norm3(r1Xr0) if nr1Xr0 / nr0 > epsilon + cross3!(r1Xr2, r1, r2) + nr1 = norm3(r1) + nr2 = norm3(r2) + @inbounds for k in 1:3 + r1r2norm[k] = r1[k]/nr1 - r2[k]/nr2 + end nr1Xr2 = norm3(r1Xr2) coeff = (gamma / (4π)) / (nr1Xr2^2) * dot3(r0, r1r2norm) @inbounds for k in 1:3 @@ -144,7 +143,6 @@ as implemented in KiteAeroDyn". v_a, work_vectors ) - r0 = work_vectors[1] r1 = work_vectors[2] r2 = work_vectors[3] r_perp = work_vectors[4] @@ -153,34 +151,29 @@ as implemented in KiteAeroDyn". r2Xr0 = work_vectors[7] normr1r2 = work_vectors[8] - r0 .= filament.x2 .- filament.x1 + r0 = filament.r0 + nr0 = filament.length r1 .= XVP .- filament.x1 r2 .= XVP .- filament.x2 - # Vector perpendicular to core radius - nr0 = norm3(r0) nr0sq = nr0 * nr0 d_r1_r0 = dot3(r1, r0) - @inbounds for k in 1:3 - r_perp[k] = d_r1_r0 * r0[k] / nr0sq - end - # Cut-off radius - epsilon = sqrt(4 * ALPHA0 * NU * norm3(r_perp) / v_a) + # Cut-off radius. The perpendicular component has length |r1.r0|/|r0|, so the + # vector itself is only needed inside the core. + epsilon = sqrt(4 * ALPHA0 * NU * abs(d_r1_r0) / nr0 / v_a) - cross3!(r1Xr2, r1, r2) cross3!(r1Xr0, r1, r0) - cross3!(r2Xr0, r2, r0) - - nr1 = norm3(r1) - nr2 = norm3(r2) - @inbounds for k in 1:3 - normr1r2[k] = r1[k]/nr1 - r2[k]/nr2 - end # Check point location relative to filament nr1Xr0 = norm3(r1Xr0) if nr1Xr0 / nr0 > epsilon + cross3!(r1Xr2, r1, r2) + nr1 = norm3(r1) + nr2 = norm3(r2) + @inbounds for k in 1:3 + normr1r2[k] = r1[k]/nr1 - r2[k]/nr2 + end nr1Xr2 = norm3(r1Xr2) coeff = (gamma / (4π)) / (nr1Xr2^2) * dot3(r0, normr1r2) @inbounds for k in 1:3 @@ -192,6 +185,7 @@ as implemented in KiteAeroDyn". # Project onto core radius — reuse r_perp, normr1r2 r1_proj = r_perp r2_proj = normr1r2 + cross3!(r2Xr0, r2, r0) nr2Xr0 = norm3(r2Xr0) d_r2_r0 = dot3(r2, r0) @inbounds for k in 1:3 @@ -264,22 +258,19 @@ function velocity_3D_trailing_vortex_semiinfinite!( work_vectors ) r1 = work_vectors[1] - r_perp = work_vectors[2] r1XVf = work_vectors[3] GAMMA = -GAMMA * filament.filament_direction r1 .= XVP .- filament.x1 - # Calculate core radius + # Core radius. `r_perp` is `(r1.Vf) Vf`, so its length is `|r1.Vf| |Vf|` and + # the vector itself is only needed inside the core. d_r1_Vf = dot3(r1, Vf) - @inbounds for k in 1:3 - r_perp[k] = d_r1_Vf * Vf[k] - end - epsilon = sqrt(4 * ALPHA0 * NU * norm3(r_perp) / v_a) + nVf = norm3(Vf) + epsilon = sqrt(4 * ALPHA0 * NU * abs(d_r1_Vf) * nVf / v_a) cross3!(r1XVf, r1, Vf) nr1XVf = norm3(r1XVf) - nVf = norm3(Vf) nr1 = norm3(r1) if nr1XVf / nVf > epsilon K = GAMMA / (4π) / (nr1XVf^2) * (1 + d_r1_Vf / nr1) diff --git a/src/solver.jl b/src/solver.jl index d008e544..3fd88574 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -130,9 +130,9 @@ Main solver structure for the Vortex Step Method.See also: [solve](@ref) (the conservative envelope from the paper) ## Additional settings -- `type_initial_gamma_distribution`::InitialGammaDistribution = ELLIPTIC: see: [InitialGammaDistribution](@ref) +- `type_initial_gamma_distribution`::InitialGammaDistribution = ZEROS: see: [InitialGammaDistribution](@ref) - `use_gamma_prev`::Bool = true: reuse provided previous gamma as initial guess when available -- `core_radius_fraction`::Float64 = 1e-20: +- `core_radius_fraction`::Float64 = 0.05: - mu::Float64 = 1.81e-5: Dynamic viscosity [N·s/m²] - `is_only_f_and_gamma_output`::Bool = false: Whether to only output f and gamma - `reference_point`::MVec3 = [0.0, 0.0, 0.0]: Moment reference point in body frame From aa57cf8974d174510ce2c441a8643967f27e2094 Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 17 Aug 2026 00:38:26 +0200 Subject: [PATCH 3/7] Add changelog entries for the filament kernel work and doc fixes Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index ba4f401c..61d74ed6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,12 +12,20 @@ them when `table_format` differs from what the directory holds, so a dataset changes format without re-running the airfoil solver that produced it. +### Fixed +- The `Solver` docstring listed `type_initial_gamma_distribution` as `ELLIPTIC` and + `core_radius_fraction` as `1e-20`; the defaults are `ZEROS` and `0.05`. + ### Changed - BREAKING: `BodyAerodynamics.AIC` is stored as `(n_panels, n_panels, 3)` instead of `(3, n_panels, n_panels)`, so each component slice `AIC[:, :, k]` is contiguous and the induced-velocity products reach BLAS `gemv` instead of the generic fallback. `solve!` is 3.1–5.3× faster (n=120, VSM: 26.1 ms → 4.9 ms inviscid, 21.1 ms → 6.0 ms with polars). Code reading `AIC[k, i, j]` must become `AIC[i, j, k]`. +- The filament induced-velocity kernels reuse the `r0`/`length` each `BoundFilament` + already stores, take the core-radius cutoff from `|r1.r0|/|r0|` without forming the + perpendicular vector, and defer the cross products that only one branch reads. Output + is bit-identical; `solve!` is a further 1.47-1.55x faster. - `read_node_table` parses into a preallocated matrix instead of `reduce(vcat, …)` over a generator, which was quadratic in the row count: ~21× faster on a 16 MB surface table (2.49 s → 0.12 s), benefiting every existing dataset. From e63eb0b2b10f37daeb6af6dcb43682f1ae72e739 Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 17 Aug 2026 09:35:05 +0200 Subject: [PATCH 4/7] Document core_radius_fraction and its Damiani reference The Solver docstring carried the SolverSettings defaults rather than its own and left core_radius_fraction undescribed. 0.05 is the upstream VSM default, attributed there to Damiani et al. (2019). Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 7 +++++-- data/pyramid_model/vsm_settings.yaml | 2 +- src/solver.jl | 3 ++- 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 61d74ed6..24a0dff6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,8 +13,11 @@ changes format without re-running the airfoil solver that produced it. ### Fixed -- The `Solver` docstring listed `type_initial_gamma_distribution` as `ELLIPTIC` and - `core_radius_fraction` as `1e-20`; the defaults are `ZEROS` and `0.05`. +- The `Solver` docstring quoted the `SolverSettings` defaults for + `type_initial_gamma_distribution` and `core_radius_fraction` (`ELLIPTIC`, `1e-20`) + instead of its own (`ZEROS`, `0.05`), and never said what `core_radius_fraction` + measures. It now documents the `Solver` defaults and cites Damiani et al. (2019) for + the 0.05 cut-off. The two structs still disagree on both values. ### Changed - BREAKING: `BodyAerodynamics.AIC` is stored as `(n_panels, n_panels, 3)` instead of diff --git a/data/pyramid_model/vsm_settings.yaml b/data/pyramid_model/vsm_settings.yaml index 74781d99..c3db7fee 100644 --- a/data/pyramid_model/vsm_settings.yaml +++ b/data/pyramid_model/vsm_settings.yaml @@ -68,7 +68,7 @@ solver_settings: artificial_damping: false # Enable artificial damping for unstable cases k2: 0.0 # 2nd-order damping coefficient k4: 0.0 # 4th-order damping coefficient - core_radius_fraction: 0.05 # Vortex core radius (fraction of chord) + core_radius_fraction: 0.05 # Vortex core radius (fraction of filament length) # --- Initial Conditions --- type_initial_gamma_distribution: ELLIPTIC # Starting circulation distribution diff --git a/src/solver.jl b/src/solver.jl index 3fd88574..1d1048cf 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -132,7 +132,8 @@ Main solver structure for the Vortex Step Method.See also: [solve](@ref) ## Additional settings - `type_initial_gamma_distribution`::InitialGammaDistribution = ZEROS: see: [InitialGammaDistribution](@ref) - `use_gamma_prev`::Bool = true: reuse provided previous gamma as initial guess when available -- `core_radius_fraction`::Float64 = 0.05: +- `core_radius_fraction`::Float64 = 0.05: Bound vortex core cut-off, as a fraction of the + filament length, following Damiani et al. (2019) - mu::Float64 = 1.81e-5: Dynamic viscosity [N·s/m²] - `is_only_f_and_gamma_output`::Bool = false: Whether to only output f and gamma - `reference_point`::MVec3 = [0.0, 0.0, 0.0]: Moment reference point in body frame From 8a6e9f305fea6efc8b2edcd9e80568469e19daba Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 17 Aug 2026 10:06:18 +0200 Subject: [PATCH 5/7] Align SolverSettings defaults with Solver core_radius_fraction 1e-20 -> 0.05 and type_initial_gamma_distribution ELLIPTIC -> ZEROS. 0.05 is the upstream awegroup/Vortex-Step-Method default, attributed there to Damiani et al. (2019), "A Vortex Step Method for Nonlinear Airfoil Polar Data as Implemented in KiteAeroDyn". At 1e-20 the Biot-Savart singularity guard never engaged. Coefficients are unchanged for well-separated geometry: the guard triggers only when a control point falls within 5% of a filament length of a bound vortex, and that ratio is 0.5*n/AR for a wing. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 10 +++++++++- data/TUDELFT_V3_KITE/vsm_settings.yaml | 2 +- data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml | 2 +- src/settings.jl | 10 +++++----- 4 files changed, 16 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 24a0dff6..1f71e5b7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,9 +17,17 @@ `type_initial_gamma_distribution` and `core_radius_fraction` (`ELLIPTIC`, `1e-20`) instead of its own (`ZEROS`, `0.05`), and never said what `core_radius_fraction` measures. It now documents the `Solver` defaults and cites Damiani et al. (2019) for - the 0.05 cut-off. The two structs still disagree on both values. + the 0.05 cut-off. ### Changed +- `SolverSettings` now defaults to the same values as `Solver`: `core_radius_fraction` + `1e-20` → `0.05` and `type_initial_gamma_distribution` `ELLIPTIC` → `ZEROS`. The 0.05 + bound vortex core cut-off follows Damiani et al. (2019), "A Vortex Step Method for + Nonlinear Airfoil Polar Data as Implemented in KiteAeroDyn", and matches the upstream + `awegroup/Vortex-Step-Method` default; at `1e-20` the Biot-Savart singularity guard + never engaged. Coefficients are unchanged for well-separated geometry, since the guard + only triggers where a control point falls within 5% of a filament length of a bound + vortex, and every settings file shipped in `data/` sets both keys explicitly. - BREAKING: `BodyAerodynamics.AIC` is stored as `(n_panels, n_panels, 3)` instead of `(3, n_panels, n_panels)`, so each component slice `AIC[:, :, k]` is contiguous and the induced-velocity products reach BLAS `gemv` instead of the generic fallback. diff --git a/data/TUDELFT_V3_KITE/vsm_settings.yaml b/data/TUDELFT_V3_KITE/vsm_settings.yaml index d5fd4c9e..86683aa2 100644 --- a/data/TUDELFT_V3_KITE/vsm_settings.yaml +++ b/data/TUDELFT_V3_KITE/vsm_settings.yaml @@ -69,7 +69,7 @@ solver_settings: artificial_damping: false # Enable artificial damping for unstable cases k2: 0.1 # 2nd-order damping coefficient k4: 0.0 # 4th-order damping coefficient - core_radius_fraction: 0.05 # Vortex core radius (fraction of chord) + core_radius_fraction: 0.05 # Vortex core radius (fraction of filament length) # --- Initial Conditions --- type_initial_gamma_distribution: ELLIPTIC # Starting circulation distribution diff --git a/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml b/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml index 4a2308cb..e29eb6d5 100644 --- a/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml +++ b/data/TUDELFT_V3_KITE/vsm_settings_coarse.yaml @@ -73,7 +73,7 @@ solver_settings: artificial_damping: false # Enable artificial damping for unstable cases k2: 0.1 # 2nd-order damping coefficient k4: 0.0 # 4th-order damping coefficient - core_radius_fraction: 1e-20 # Vortex core radius (fraction of chord) + core_radius_fraction: 1e-20 # Vortex core radius (fraction of filament length) # --- Initial Conditions --- type_initial_gamma_distribution: ELLIPTIC # Starting circulation distribution diff --git a/src/settings.jl b/src/settings.jl index 83b05c90..33fa36d3 100644 --- a/src/settings.jl +++ b/src/settings.jl @@ -79,11 +79,11 @@ Solver configuration, used within [`VSMSettings`](@ref). (default `0.035`) - `type_initial_gamma_distribution`: [`ELLIPTIC`](@ref InitialGammaDistribution) or `ZEROS` - (default `ELLIPTIC`) + (default `ZEROS`) - `use_gamma_prev`: Reuse provided previous gamma as initial guess when available (default `true`) -- `core_radius_fraction`: Vortex core radius fraction - (default `1e-20`) +- `core_radius_fraction`: Bound vortex core cut-off, as a fraction of the filament + length, following Damiani et al. (2019) (default `0.05`) - `mu`: Dynamic viscosity (N*s/m^2) (default `1.81e-5`) - `calc_only_f_and_gamma`: Only output forces and circulation (default `false`) @@ -104,9 +104,9 @@ Solver configuration, used within [`VSMSettings`](@ref). k4::Float64 = 0.0 # artificial damping parameter is_with_artificial_viscosity::Bool = false # Li/Gaunaa post-stall artificial viscosity artificial_viscosity_factor::Float64 = 0.035 # viscosity scaling coefficient k - type_initial_gamma_distribution::InitialGammaDistribution = ELLIPTIC # see: [InitialGammaDistribution](@ref) + type_initial_gamma_distribution::InitialGammaDistribution = ZEROS # see: [InitialGammaDistribution](@ref) use_gamma_prev::Bool = true # if false, always reinitialize gamma from type_initial_gamma_distribution - core_radius_fraction::Float64 = 1e-20 + core_radius_fraction::Float64 = 0.05 mu::Float64 = 1.81e-5 # dynamic viscosity [N·s/m²] calc_only_f_and_gamma::Bool=false # whether to only output f and gamma correct_aoa::Bool=false # perform aoa correction From ae7294d1fdb1169c49723f0eff8eaa6a825fee3e Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 17 Aug 2026 10:41:23 +0200 Subject: [PATCH 6/7] Keep test-env precompilation lazy on macOS and Windows Julia 1.12.7 precompiles the whole test environment before running any test, so MakieControlPlots and GLMakie load ahead of the Sys.islinux() guard in runtests.jl. Only Linux sets up xvfb, so GLMakie fails to open a GL context there: "GLFWError (FORMAT_UNAVAILABLE): NSGL: Failed to find a suitable pixel format". Same runner image, same GLMakie/GLFW versions as the last green run on 1.12.6. Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/CI.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 99f83366..0b6f5dae 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -64,6 +64,10 @@ jobs: - uses: julia-actions/julia-runtest@v1 env: BUILD_IS_PRODUCTION_BUILD: ${{ matrix.build_is_production_build }} + # Only Linux gets a virtual display, so runtests.jl skips the plotting tests + # elsewhere. Julia 1.12.7 precompiles the whole test environment up front, + # which loads GLMakie before that guard runs; keep it lazy off Linux. + JULIA_PKG_PRECOMPILE_AUTO: ${{ runner.os == 'Linux' && '1' || '0' }} with: coverage: false prefix: ${{ runner.os == 'Linux' && 'xvfb-run -a' || '' }} From 2d6ed041fce9c2fb305591597c736e21faf739a4 Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 17 Aug 2026 11:35:11 +0200 Subject: [PATCH 7/7] Drop the BREAKING prefix on the AIC layout entry AIC is an internal buffer; its only accessor, calculate_AIC_matrices!, is documented under private functions. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1f71e5b7..2e73b25a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,7 +28,7 @@ never engaged. Coefficients are unchanged for well-separated geometry, since the guard only triggers where a control point falls within 5% of a filament length of a bound vortex, and every settings file shipped in `data/` sets both keys explicitly. -- BREAKING: `BodyAerodynamics.AIC` is stored as `(n_panels, n_panels, 3)` instead of +- `BodyAerodynamics.AIC` is stored as `(n_panels, n_panels, 3)` instead of `(3, n_panels, n_panels)`, so each component slice `AIC[:, :, k]` is contiguous and the induced-velocity products reach BLAS `gemv` instead of the generic fallback. `solve!` is 3.1–5.3× faster (n=120, VSM: 26.1 ms → 4.9 ms inviscid, 21.1 ms →