Truncate the SVD in the least-squares Kulfan fit, and warn when it does, so crowded stations give bounded weights

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.

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

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]

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(