diff --git a/CHANGELOG.md b/CHANGELOG.md index 3c2babbf..13e1e426 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,10 @@ - The `VSMSolution` docstring gives `lift_dist`, `drag_dist` and `panel_moment_dist` in the per-unit-span units they hold, [N/m] and [Nm/m], instead of [N] and [Nm]. +- `fit_kulfan_parameters` with `LeastSquaresFit` drops singular values below `1e-4` + times the largest and warns when it does, so a contour whose stations crowd into a + narrow band of the chord gets bounded weights instead of ones that resample it to 1e4 + scale. Fits of well-spread stations are unchanged. - Inside its vortex core, `velocity_3D_trailing_vortex!` induces an azimuthal velocity instead of a radial one. Only points within the millimetre-scale Oseen core of a panel's chordwise trailing segment were affected. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 5f7eec56..5f954cb1 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -150,6 +150,7 @@ chord_residual! bernstein_basis class_function leading_edge_basis +truncated_least_squares compute_optimal_x_points normalize_airfoil get_lower_upper diff --git a/src/airfoil_aero/kulfan.jl b/src/airfoil_aero/kulfan.jl index dce90044..7d7f2c30 100644 --- a/src/airfoil_aero/kulfan.jl +++ b/src/airfoil_aero/kulfan.jl @@ -106,6 +106,21 @@ function leading_edge_basis(x::AbstractVector{T}, n_weights::Int) where T return x .* max.(1 .- x, zero(T)).^(n_weights + 0.5) end +const KULFAN_FIT_RTOL = 1e-4 # [-] singular values dropped below this times the largest + +""" + truncated_least_squares(A, b) + +Minimum-norm least-squares solution of `A * coeffs = b` with the singular values of `A` +below `KULFAN_FIT_RTOL` times the largest dropped. Returns `(coeffs, n_dropped)`. +""" +function truncated_least_squares(A::AbstractMatrix, b::AbstractVector) + F = svd(A) + kept = F.S .> KULFAN_FIT_RTOL * F.S[1] + coeffs = F.V[:, kept] * ((F.U[:, kept]' * b) ./ F.S[kept]) + return coeffs, count(!, kept) +end + """ fit_kulfan_parameters(x::Vector, y::Vector, method::KulfanFitMethod) fit_kulfan_parameters(x::Vector, y::Vector; n_weights=8) @@ -128,7 +143,8 @@ end Least-squares fit matching AeroSandbox's `get_kulfan_parameters`: both surfaces share a single least-squares system with a shared leading-edge weight and a -trailing-edge thickness. +trailing-edge thickness. Singular values below `1e-4` times the largest are dropped, so +stations crowded into part of the chord give bounded weights, and a warning says so. """ function fit_kulfan_parameters(x::Vector{T}, y::Vector{T}, method::LeastSquaresFit) where T @@ -144,14 +160,16 @@ function fit_kulfan_parameters(x::Vector{T}, y::Vector{T}, te_col = ifelse.(is_upper, xv ./ 2, .-xv ./ 2) A = hcat((.!is_upper) .* CS, is_upper .* CS, le_col, te_col) - coeffs = A \ y_norm + coeffs, n_dropped = truncated_least_squares(A, y_norm) TE_thickness = coeffs[end] if TE_thickness < 0 - A = hcat((.!is_upper) .* CS, is_upper .* CS, le_col) - coeffs = A \ y_norm + coeffs, n_dropped = truncated_least_squares(A[:, 1:end-1], y_norm) TE_thickness = zero(T) end + n_dropped > 0 && @warn "Kulfan fit dropped $n_dropped of $(length(coeffs)) singular \ + values: the stations leave part of the shape unconstrained, as when they crowd \ + into a narrow band of the chord." lower_weights = coeffs[1:n_weights] upper_weights = coeffs[n_weights+1:2n_weights] diff --git a/test/airfoil_aero/test_airfoil_aero.jl b/test/airfoil_aero/test_airfoil_aero.jl index f306968a..65329629 100644 --- a/test/airfoil_aero/test_airfoil_aero.jl +++ b/test/airfoil_aero/test_airfoil_aero.jl @@ -43,6 +43,23 @@ end @test params.TE_thickness ≈ 0.0 atol = 1e-12 end + @testset "Fit to stations crowded into a narrow band stays airfoil-sized" begin + spread = (1 .- cos.(range(0, pi, 121))) ./ 2 + crowded = vcat(0.0, range(0.30, 0.31, 118), 1.0) + surface(weights, xs) = class_function(xs) .* (bernstein_basis(xs, 7) * weights) + y_upper = surface(fill(0.2, 8), crowded) .+ 0.005 .* sin.(40pi .* crowded) + y_lower = surface(fill(-0.1, 8), spread) + x = vcat(reverse(crowded), spread[2:end]) + y = vcat(reverse(y_upper), y_lower[2:end]) + crowded_fit = @test_logs (:warn, r"dropped \d+ of \d+ singular values") ( + fit_kulfan_parameters(x, y)) + _, y_fit = kulfan_to_coordinates(crowded_fit) + @test maximum(abs, y_fit) < 2 * maximum(abs, y) + + y_spread = vcat(reverse(surface(fill(0.2, 8), spread)), y_lower[2:end]) + @test_logs fit_kulfan_parameters(vcat(reverse(spread), spread[2:end]), y_spread) + end + @testset "Shrink-wrap encloses points with clearance" begin xn, yn, _ = normalize_airfoil(collect(float.(xr)), collect(float.(yr))) cloud_to_wrap(xw, yw) = minimum(