From 78b0681842ce80972a5ad26e9279ec1f261a18e2 Mon Sep 17 00:00:00 2001 From: Bart Date: Tue, 25 Aug 2026 12:38:32 +0200 Subject: [PATCH 1/8] Regenerate polars from the deformed shape instead of a flap angle MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A trailing edge carried by a chordwise beam bends into a smooth camber shape. Projecting that onto one hinge angle throws most of it away, and which angle you get depends on which node pair you measure between. Add the pieces to skip the projection. `TAYLOR` is a per-panel polynomial of arbitrary order in `α - α_ref`, the local expansion a live source refits each solve; it is valid only inside its window, so it ignores δ and is skipped by the stall-angle scan. `KulfanBasis` and `deform_kulfan` deform a fixed Kulfan fit by a matvec against a constant CST basis — never by refitting, which is non-unique enough that the same shape fitted twice moves Cl by more than the deformation does. `control_point_deflection` resamples a deflection given at arbitrary chord fractions onto that basis, so beam nodes today and membrane nodes later feed the same path. `LivePolars` then evaluates NeuralFoil a few angles either side of each panel's own α, in one batched forward pass, and least-squares fits the coefficients. Order and window belong together: order 2 over ±4° tracks a full-polar solve to a few tenths of a percent and bends with the stall knee, while higher orders diverge as soon as the solve steps outside. `polar_drift` reports how far it has, so the caller can refit before trusting it. Costs 2.8 ms for 20 panels, 7.1 ms for 80. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 17 +++ docs/src/functions.md | 5 + docs/src/private_functions.md | 4 + docs/src/private_types.md | 3 + src/VortexStepMethod.jl | 17 ++- src/airfoil_aero/AirfoilAero.jl | 7 +- src/airfoil_aero/deform.jl | 92 ++++++++++++++++ src/airfoil_aero/live_polar.jl | 158 +++++++++++++++++++++++++++ src/airfoil_aero/neuralfoil.jl | 106 ++++++++++++------ src/body_aerodynamics.jl | 7 ++ src/panel.jl | 41 +++++++ src/wing_geometry.jl | 14 +++ test/airfoil_aero/test_live_polar.jl | 111 +++++++++++++++++++ test/runtests.jl | 1 + 14 files changed, 542 insertions(+), 41 deletions(-) create mode 100644 src/airfoil_aero/deform.jl create mode 100644 src/airfoil_aero/live_polar.jl create mode 100644 test/airfoil_aero/test_live_polar.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index 038f8be6..9dac1907 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,23 @@ ## VortexStepMethod v4.2.0 2026-08-25 ### Added +- `TAYLOR` aero model: a per-panel polynomial of arbitrary order in `α - α_ref` + [rad], the local expansion a live polar source refits every solve. Sections + carry `aero_data = (alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs)` and blend + spanwise like any other model; `set_taylor_polar!` rewrites a panel's fit in + place. It is valid only inside its fit window, so it ignores `delta` and is + skipped by the stall-angle scan. +- Live in-memory polars in `AirfoilAero`: `KulfanBasis` and `deform_kulfan` + deform a fixed Kulfan fit analytically (a matvec against a constant CST basis, + never a refit, which is non-unique); `control_point_deflection` resamples a + deflection given at arbitrary chord fractions — beam nodes, membrane nodes — + onto that basis. `LivePolars` and `refresh_live_polars!` then evaluate + NeuralFoil at a few angles per panel in one batched forward pass and + least-squares fit each panel's `TAYLOR` polar, so a chordwise deformation + reaches the aerodynamics as a shape change rather than a flap angle. + `polar_drift` reports how far the solve has left the fit window. +- `prepare_inputs` accepts one Kulfan shape per case, so a whole wing goes + through NeuralFoil in a single forward pass. - Shared panel aerodynamics (`src/panel_aerodynamics.jl`): the per-panel physics written once as pure, branch-free, number-type-generic functions of the section geometry and the flow — `panel_axes`, `panel_inflow`, `panel_force_directions`, diff --git a/docs/src/functions.md b/docs/src/functions.md index 6562b163..380b8523 100644 --- a/docs/src/functions.md +++ b/docs/src/functions.md @@ -41,6 +41,11 @@ deform_section analyze_section analyze_sweep neuralfoil_aero +deform_kulfan +control_point_deflection +panel_kulfan_parameters +refresh_live_polars! +polar_drift generate_aero_matrices generate_polar_from_coordinates generate_polar_from_dat diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index b48ed3e8..3453a80e 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -98,6 +98,7 @@ update_non_deformed_sections! ### Aerodynamic data and Cp ```@docs calculate_new_aero_data +set_taylor_polar! assemble_polar_matrix load_matrix_polar_data read_aero_matrix @@ -153,8 +154,11 @@ smooth_turning! load_neuralfoil_model neuralfoil_section neuralfoil_fused_output +fused_output +decode_coefficients nn_forward prepare_inputs +fill_case_input! flip_inputs flip_outputs squared_mahalanobis_distance diff --git a/docs/src/private_types.md b/docs/src/private_types.md index 5fca5ebe..21e223d3 100644 --- a/docs/src/private_types.md +++ b/docs/src/private_types.md @@ -24,6 +24,9 @@ CurrentModule = VortexStepMethod.AirfoilAero ``` ```@docs KulfanParameters +KulfanBasis +LivePolarSettings +LivePolars NeuralFoilModel NeuralFoilResult ``` diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 087a346a..93661842 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -36,7 +36,7 @@ export calculate_projected_area, calculate_span export MVec3 export LLT, Model, VSM -export AeroModel, INVISCID, POLY, LEI_AIRFOIL_BREUKELS, POLAR_MATRICES, POLAR_VECTORS +export AeroModel, INVISCID, POLY, LEI_AIRFOIL_BREUKELS, POLAR_MATRICES, POLAR_VECTORS, TAYLOR export BILLOWING, COSINE, LINEAR, PanelDistribution, SPLIT_PROVIDED, UNCHANGED export ELLIPTIC, InitialGammaDistribution, ZEROS export FAILURE, FEASIBLE, INFEASIBLE, SolverStatus @@ -237,7 +237,7 @@ Enumeration of the implemented wing types. @enum WingType RECTANGULAR CURVED ELLIPTICAL """ - AeroModel `POLY` `POLAR_VECTORS` `POLAR_MATRICES` `INVISCID` + AeroModel `POLY` `POLAR_VECTORS` `POLAR_MATRICES` `INVISCID` `TAYLOR` Enumeration of the implemented aerodynamic models. See also: [AeroData](@ref) @@ -247,6 +247,10 @@ Enumeration of the implemented aerodynamic models. See also: [AeroData](@ref) - `POLAR_VECTORS`: Polar vectors as function of alpha (lookup tables with interpolation) - `POLAR_MATRICES`: Polar matrices as function of alpha and delta (lookup tables with interpolation) - INVISCID +- `TAYLOR`: order-N polynomial in `α - α_ref` [rad] per panel, the local expansion a + live polar source refits each solve, see [`refresh_live_polars!`](@ref + VortexStepMethod.AirfoilAero.refresh_live_polars!). Valid only inside its fit + window, so it carries no stall model and ignores `delta`. `LEI_AIRFOIL_BREUKELS` is a deprecated alias of `POLY`. @@ -257,6 +261,7 @@ where `alpha` is the angle of attack, `delta` is trailing edge angle. POLAR_VECTORS POLAR_MATRICES INVISCID + TAYLOR end """ @@ -330,7 +335,8 @@ abstract type AbstractWing{T} end Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}}, - Tuple{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}} + Tuple{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, + Tuple{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}} } Union of different definitions of the aerodynamic properties of a wing section. See also: [AeroModel](@ref) @@ -338,6 +344,8 @@ Union of different definitions of the aerodynamic properties of a wing section. - (`cl_coeffs`, `cd_coeffs`, `cm_coeffs`) α-polynomial coefficients for `POLY` - (`alpha_range`, `cl_vector`, `cd_vector`, `cm_vector`) for `POLAR_VECTORS` - (`alpha_range`, `delta_range`, `cl_matrix`, `cd_matrix`, `cm_matrix`) for `POLAR_MATRICES` + - (`alpha_ref`, `cl_coeffs`, `cd_coeffs`, `cm_coeffs`) for `TAYLOR`, the expansion + point [rad] and the ascending coefficients of the polynomial in `α - alpha_ref` where `alpha` is the angle of attack [rad], `delta` is trailing edge angle [rad], `cl` the lift coefficient, `cd` the drag coefficient and `cm` the pitching moment coefficient. The camber of a kite refers to @@ -349,7 +357,8 @@ const AeroData = Union{ Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}}, - Tuple{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}} + Tuple{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, + Tuple{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}} } const PACKAGE_ROOT = normpath(joinpath(@__DIR__, "..")) diff --git a/src/airfoil_aero/AirfoilAero.jl b/src/airfoil_aero/AirfoilAero.jl index 531a5ad2..452f775a 100644 --- a/src/airfoil_aero/AirfoilAero.jl +++ b/src/airfoil_aero/AirfoilAero.jl @@ -8,11 +8,13 @@ using NPZ using Xfoil using Printf: @sprintf using ..VortexStepMethod: SectionAero, interpolate_matrix_nans!, delta_suffix, - write_node_rows + write_node_rows, section_surface, set_taylor_polar! include("kulfan.jl") +include("deform.jl") include("shrink_wrap.jl") include("neuralfoil.jl") +include("live_polar.jl") include("poly.jl") include("airfoil_solvers/common.jl") include("airfoil_solvers/xfoil_solver.jl") @@ -29,6 +31,9 @@ export ShrinkWrap, shrink_wrap export fit_kulfan_parameters, kulfan_to_coordinates export NeuralFoilModel, NeuralFoilResult, load_neuralfoil_model export neuralfoil_aero, neuralfoil_section +export KulfanBasis, deform_kulfan, control_point_deflection +export LivePolarSettings, LivePolars, panel_kulfan_parameters +export refresh_live_polars!, polar_drift export AbstractAirfoilSolver, XFoilSolver, NeuralFoilSolver export SectionSolution, DeformedSection, deform_section, analyze_section, analyze_sweep export create_2d_polars, generate_aero_matrices, generate_section_aero, lei_poly_coeffs diff --git a/src/airfoil_aero/deform.jl b/src/airfoil_aero/deform.jl new file mode 100644 index 00000000..52a165f5 --- /dev/null +++ b/src/airfoil_aero/deform.jl @@ -0,0 +1,92 @@ +""" + KulfanBasis(; n_stations=60, n_weights=8) + +The fixed CST basis a shape deformation is projected onto: chord stations `x` in +`[0, 1]` and the pseudoinverse `projection` of the class-times-Bernstein matrix +`C(x)·B(x)` those stations span. + +CST is linear in its weights, so a surface displacement is a matvec against this +constant matrix — never a refit. Refitting inside a loop is not an option: the Kulfan +fit is non-unique, so the same shape fitted twice returns weight vectors differing by +more than the deformation signal, which reaches the polars as frame-to-frame jitter. + +`C(1) = 0`, so the parameterisation cannot put the trailing edge off the chord line. +Deflections must therefore be measured **against the deformed chord**, with the chord +rotation and stretch taken from the leading- and trailing-edge points; what is left is +the representable residual. +""" +struct KulfanBasis + "Chord stations the deflection is sampled on, ascending in `[0, 1]`." + x::Vector{Float64} + "`n_weights × length(x)` pseudoinverse of `C(x)·B(x)`, mapping deflection to weights." + projection::Matrix{Float64} + "Number of CST weights per surface, matching the airfoil being deformed." + n_weights::Int +end + +function KulfanBasis(; n_stations::Int=60, n_weights::Int=8) + n_stations > n_weights || throw(ArgumentError( + "KulfanBasis needs more stations than weights; got $n_stations and $n_weights.")) + theta = range(0, pi, n_stations) + x = @. (1 - cos(theta)) / 2 + shape = class_function(x) .* bernstein_basis(x, n_weights - 1) + return KulfanBasis(collect(x), pinv(shape), n_weights) +end + +""" + deform_kulfan(basis, base, upper_deflection, lower_deflection) -> KulfanParameters + deform_kulfan(basis, base, camber) -> KulfanParameters + +Add a surface deflection to a fixed set of Kulfan parameters, both deflections sampled +on `basis.x` and normalized by chord. The three-argument form applies one camber +deflection to both surfaces, leaving the thickness distribution untouched; the +four-argument form deforms the surfaces independently, which is what a double-skin +membrane with its own upper and lower control points needs. + +The leading-edge weight and the trailing-edge thickness are carried over unchanged: +they are the two shape freedoms a chord-referenced deflection cannot resolve. +""" +function deform_kulfan(basis::KulfanBasis, base::KulfanParameters, + upper_deflection::AbstractVector, + lower_deflection::AbstractVector) + length(base.upper_weights) == basis.n_weights || throw(ArgumentError( + "KulfanBasis has $(basis.n_weights) weights, airfoil has " * + "$(length(base.upper_weights)).")) + (length(upper_deflection) == length(basis.x) && + length(lower_deflection) == length(basis.x)) || throw(ArgumentError( + "Deflections must be sampled on the basis' $(length(basis.x)) stations.")) + return KulfanParameters(base.upper_weights .+ basis.projection * upper_deflection, + base.lower_weights .+ basis.projection * lower_deflection, + base.leading_edge_weight, base.TE_thickness) +end + +deform_kulfan(basis::KulfanBasis, base::KulfanParameters, + camber::AbstractVector) = deform_kulfan(basis, base, camber, camber) + +""" + control_point_deflection(basis, fractions, deflections) -> Vector{Float64} + +Resample a deflection given at arbitrary chord `fractions` onto `basis.x`, by linear +interpolation, held flat outside the sampled range. This is the generic bridge from +control points — beam nodes now, membrane nodes later — to the CST basis: the caller +only has to say where along the chord each point sits and how far off the chord line it +moved, both normalized by the local chord. + +`fractions` need not be sorted, but must not repeat a station. +""" +function control_point_deflection(basis::KulfanBasis, fractions::AbstractVector, + deflections::AbstractVector) + length(fractions) == length(deflections) || throw(ArgumentError( + "control_point_deflection: $(length(fractions)) fractions for " * + "$(length(deflections)) deflections.")) + order = sortperm(collect(float.(fractions))) + knots = collect(float.(fractions))[order] + values = collect(float.(deflections))[order] + if length(knots) < 2 + return fill(isempty(values) ? 0.0 : values[1], length(basis.x)) + end + allunique(knots) || throw(ArgumentError( + "control_point_deflection: repeated chord fraction in $knots.")) + interp = linear_interpolation(knots, values; extrapolation_bc=Flat()) + return [interp(xi) for xi in basis.x] +end diff --git a/src/airfoil_aero/live_polar.jl b/src/airfoil_aero/live_polar.jl new file mode 100644 index 00000000..d1411c8c --- /dev/null +++ b/src/airfoil_aero/live_polar.jl @@ -0,0 +1,158 @@ +""" + LivePolarSettings(; order=2, half_window=deg2rad(4), n_samples=5, + model_size="xlarge", weights_dir=nothing, n_crit=9.0) + +How a live polar is sampled and fitted. Every solve, each panel's deformed shape is +evaluated at `n_samples` angles of attack spanning `α_ref ± half_window` and a +polynomial of `order` is least-squares fitted through them, becoming the panel's +`TAYLOR` polar. + +Order and window belong together. Order 2 over `±4°` tracks a full-polar solve to +within a few tenths of a percent and bends with the stall knee where a straight line +cuts across it; higher orders carry large coefficients that diverge as soon as the +solve steps outside the window, which is why the window is a correctness limit and not +an accuracy knob — see [`polar_drift`](@ref). + +Fitting is by least squares over the window, never by finite differences at a tabulated +source's own knot spacing: a piecewise-linear polar has zero curvature inside a segment +and a spike at every knot, so a difference at the knot spacing returns pure +discretization artefact. +""" +@with_kw struct LivePolarSettings + "Polynomial order in `α - α_ref`; 2 unless the window is narrowed to match." + order::Int = 2 + "Half width [rad] of the fit window, also the drift the solve may not leave." + half_window::Float64 = deg2rad(4.0) + "Angles of attack sampled per panel per refresh." + n_samples::Int = 5 + "NeuralFoil network size." + model_size::String = "xlarge" + "Directory holding the network weights; `nothing` takes the packaged ones." + weights_dir::Union{Nothing, String} = nothing + "Critical amplification factor of the transition model." + n_crit::Float64 = 9.0 +end + +""" + LivePolars(base; settings=LivePolarSettings(), n_stations=60) + +Live polar source for a wing whose panels each carry one undeformed airfoil in `base`. +Holds the fixed CST basis, the per-panel expansion points and the scratch the refresh +reuses, so a steady-state refresh allocates nothing beyond the network's own forward +pass. Drive it with [`refresh_live_polars!`](@ref). +""" +mutable struct LivePolars + "Sampling and fitting configuration." + settings::LivePolarSettings + "The fixed CST basis every panel's deflection is projected onto." + basis::KulfanBasis + "Undeformed Kulfan parameters per panel." + base::Vector{KulfanParameters} + "Deformed Kulfan parameters per panel, rewritten every refresh." + deformed::Vector{KulfanParameters} + "Expansion point [rad] the last fit was built about, per panel." + alpha_ref::Vector{Float64} + "Sample offsets [rad] off the expansion point, shared by every panel." + offsets::Vector{Float64} + "`(order + 1) × n_samples` least-squares solver for the shared offset grid." + fit::Matrix{Float64} + "`25 × (n_panels · n_samples)` network input scratch." + inputs::Matrix{Float32} +end + +function LivePolars(base::AbstractVector{KulfanParameters}; + settings::LivePolarSettings=LivePolarSettings(), + n_stations::Int=60) + settings.n_samples > settings.order || throw(ArgumentError( + "LivePolars needs more samples than the fit order; got " * + "$(settings.n_samples) and $(settings.order).")) + offsets = collect(range(-settings.half_window, settings.half_window, + settings.n_samples)) + vandermonde = [d^(k - 1) for d in offsets, k in 1:(settings.order + 1)] + n_panels = length(base) + return LivePolars(settings, KulfanBasis(; n_stations, + n_weights=length(base[1].upper_weights)), + collect(base), collect(base), zeros(n_panels), offsets, + pinv(vandermonde), + zeros(Float32, 25, n_panels * settings.n_samples)) +end + +""" + panel_kulfan_parameters(panels; delta=0.0) -> Vector{KulfanParameters} + +Fit the undeformed Kulfan parameters of every panel from the surface contour its +`section_aero` carries, the starting point [`LivePolars`](@ref) deforms. Fitted once at +build time, never inside the loop. +""" +function panel_kulfan_parameters(panels; delta=0.0) + return map(panels) do panel + panel.section_aero === nothing && throw(ArgumentError( + "Live polars need a section_aero contour on every panel.")) + x, y, _, _ = section_surface(panel.section_aero, 0.0, delta) + fit_kulfan_parameters(collect(float.(x)), collect(float.(y))) + end +end + +""" + polar_drift(live, alpha) -> Float64 + +How far the largest panel angle of attack has drifted off the expansion point its +polar was fitted about, as a fraction of the fit window. Above `1` a panel is being +evaluated outside the window and the fit must be rebuilt before its answer is used, so +this is the guard [`refresh_live_polars!`](@ref) leaves to the caller's solve loop. +""" +function polar_drift(live::LivePolars, alpha::AbstractVector) + length(alpha) == length(live.alpha_ref) || throw(ArgumentError( + "polar_drift: $(length(alpha)) angles for $(length(live.alpha_ref)) panels.")) + return maximum(abs.(alpha .- live.alpha_ref)) / live.settings.half_window +end + +""" + refresh_live_polars!(live, panels, alpha_ref, reynolds; deflection=nothing) + -> Float64 + +Refit every panel's polar from its current shape and write it in as a `TAYLOR` polar +(see [`set_taylor_polar!`](@ref VortexStepMethod.set_taylor_polar!)). Per panel: +deform the base airfoil by `deflection` (a chord-normalized deflection on +`live.basis.x`, or `nothing` to keep the base shape), evaluate NeuralFoil at +`alpha_ref ± half_window`, and least-squares fit the polynomial. + +`alpha_ref` [rad] and `reynolds` are per panel or scalar. Every panel and sample goes +through the network in one forward pass. Returns the lowest analysis confidence over +the batch, a value near zero meaning the deformed shape has left the region the network +was trained on. + +Call it after the mesh has been rebuilt from the structure and before the solve: a mesh +rebuild re-seeds each panel's aero from its section, which would drop the fit. +""" +function refresh_live_polars!(live::LivePolars, panels, alpha_ref, reynolds; + deflection=nothing) + n_panels = length(live.base) + length(panels) == n_panels || throw(ArgumentError( + "refresh_live_polars!: $(length(panels)) panels for $n_panels base airfoils.")) + per_panel(v) = v isa Number ? fill(float(v), n_panels) : collect(float.(v)) + alpha_vec, re_vec = per_panel(alpha_ref), per_panel(reynolds) + live.alpha_ref .= alpha_vec + + n_samples = live.settings.n_samples + for i in 1:n_panels + live.deformed[i] = isnothing(deflection) ? live.base[i] : + deform_kulfan(live.basis, live.base[i], deflection[i]) + for k in 1:n_samples + fill_case_input!(live.inputs, (i - 1) * n_samples + k, live.deformed[i], + rad2deg(alpha_vec[i] + live.offsets[k]), re_vec[i], + live.settings.n_crit, 1.0, 1.0) + end + end + + model = load_neuralfoil_model(live.settings.model_size; + weights_dir=live.settings.weights_dir) + cl, cd, cm, confidence = decode_coefficients(fused_output(live.inputs, model)) + + for i in 1:n_panels + window = ((i - 1) * n_samples + 1):(i * n_samples) + set_taylor_polar!(panels[i], alpha_vec[i], live.fit * cl[window], + live.fit * cd[window], live.fit * cm[window]) + end + return minimum(confidence) +end diff --git a/src/airfoil_aero/neuralfoil.jl b/src/airfoil_aero/neuralfoil.jl index aea9d062..ea028f57 100644 --- a/src/airfoil_aero/neuralfoil.jl +++ b/src/airfoil_aero/neuralfoil.jl @@ -165,11 +165,40 @@ function squared_mahalanobis_distance(x::AbstractMatrix, model::NeuralFoilModel) return result end +""" + fill_case_input!(x, case, params::KulfanParameters, alpha_deg, Re, n_crit, + xtr_upper, xtr_lower) + +Write one NeuralFoil input column: rows 1–18 the Kulfan shape, rows 19–25 the flow +condition (`alpha_deg` in degrees). The single point where the network's input layout +is defined, shared by every `prepare_inputs` method. +""" +function fill_case_input!(x::AbstractMatrix, case::Int, params::KulfanParameters, + alpha_deg, Re, n_crit, xtr_upper, xtr_lower) + for i in 1:8 + x[i, case] = params.upper_weights[i] + x[8 + i, case] = params.lower_weights[i] + end + x[17, case] = params.leading_edge_weight + x[18, case] = params.TE_thickness * 50 # Scale factor from NeuralFoil + x[19, case] = sind(2 * alpha_deg) + x[20, case] = cosd(alpha_deg) + x[21, case] = 1 - cosd(alpha_deg)^2 + x[22, case] = (log(Re) - 12.5) / 3.5 + x[23, case] = (n_crit - 9) / 4.5 + x[24, case] = xtr_upper + x[25, case] = xtr_lower + return nothing +end + """ prepare_inputs(params::KulfanParameters, alpha, Re; n_crit=9.0, xtr_upper=1.0, xtr_lower=1.0) + prepare_inputs(params::AbstractVector{KulfanParameters}, alpha, Re; kwargs...) -Prepare neural network inputs from Kulfan parameters and flow conditions. +Prepare neural network inputs from Kulfan parameters and flow conditions. The vector +form takes one shape per case, so a whole wing's panels go through the network in a +single forward pass. # Arguments - `params`: Kulfan CST parameters @@ -184,39 +213,27 @@ Prepare neural network inputs from Kulfan parameters and flow conditions. """ function prepare_inputs(params::KulfanParameters, alpha, Re; n_crit=9.0, xtr_upper=1.0, xtr_lower=1.0) - # Ensure vectors alpha_vec = alpha isa Number ? [alpha] : collect(alpha) - Re_vec = Re isa Number ? fill(Re, length(alpha_vec)) : collect(Re) n_cases = length(alpha_vec) + return prepare_inputs(fill(params, n_cases), alpha_vec, Re; + n_crit, xtr_upper, xtr_lower) +end - # Broadcast other parameters - n_crit_vec = n_crit isa Number ? fill(n_crit, n_cases) : collect(n_crit) - xtr_upper_vec = xtr_upper isa Number ? fill(xtr_upper, n_cases) : collect(xtr_upper) - xtr_lower_vec = xtr_lower isa Number ? fill(xtr_lower, n_cases) : collect(xtr_lower) +function prepare_inputs(params::AbstractVector{KulfanParameters}, alpha, Re; + n_crit=9.0, xtr_upper=1.0, xtr_lower=1.0) + alpha_vec = alpha isa Number ? fill(alpha, length(params)) : collect(alpha) + n_cases = length(alpha_vec) + length(params) == n_cases || throw(ArgumentError( + "prepare_inputs: $(length(params)) shapes for $n_cases angles of attack.")) + per_case(v) = v isa Number ? fill(v, n_cases) : collect(v) + Re_vec, n_crit_vec = per_case(Re), per_case(n_crit) + xtr_upper_vec, xtr_lower_vec = per_case(xtr_upper), per_case(xtr_lower) - # Build input matrix (25 x n_cases) x = zeros(Float32, 25, n_cases) - - # Kulfan parameters (same for all cases) - for i in 1:8 - x[i, :] .= params.upper_weights[i] - x[8+i, :] .= params.lower_weights[i] - end - x[17, :] .= params.leading_edge_weight - x[18, :] .= params.TE_thickness * 50 # Scale factor from NeuralFoil - - # Flow conditions (vary per case) for i in 1:n_cases - a = alpha_vec[i] - x[19, i] = sind(2 * a) - x[20, i] = cosd(a) - x[21, i] = 1 - cosd(a)^2 - x[22, i] = (log(Re_vec[i]) - 12.5) / 3.5 - x[23, i] = (n_crit_vec[i] - 9) / 4.5 - x[24, i] = xtr_upper_vec[i] - x[25, i] = xtr_lower_vec[i] + fill_case_input!(x, i, params[i], alpha_vec[i], Re_vec[i], n_crit_vec[i], + xtr_upper_vec[i], xtr_lower_vec[i]) end - return x end @@ -323,15 +340,22 @@ function neuralfoil_aero(params::KulfanParameters, alpha, Re; n_crit=9.0, xtr_upper=1.0, xtr_lower=1.0) y = neuralfoil_fused_output(params, alpha, Re; model_size, weights_dir, n_crit, xtr_upper, xtr_lower) - analysis_confidence = sigmoid.(y[1, :]) - CL = y[2, :] ./ 2 - CD = clamp.(exp.((y[3, :] .- 2) .* 2), 0.0, 1.0) - CM = y[4, :] ./ 20 - + cl, cd, cm, confidence = decode_coefficients(y) alpha_vec = alpha isa Number ? [Float64(alpha)] : Float64.(collect(alpha)) + return NeuralFoilResult(alpha_vec, cl, cd, cm, confidence) +end + +""" + decode_coefficients(y) -> (cl, cd, cm, confidence) - return NeuralFoilResult(alpha_vec, Vector{Float64}(CL), Vector{Float64}(CD), - Vector{Float64}(CM), Vector{Float64}(analysis_confidence)) +Turn a fused network output matrix into the integrated coefficients per case, undoing +NeuralFoil's output scaling. The single place that scaling is written down. +""" +function decode_coefficients(y::AbstractMatrix) + return (Vector{Float64}(y[2, :] ./ 2), + Vector{Float64}(clamp.(exp.((y[3, :] .- 2) .* 2), 0.0, 1.0)), + Vector{Float64}(y[4, :] ./ 20), + Vector{Float64}(sigmoid.(y[1, :]))) end """ @@ -343,11 +367,21 @@ top/bottom-flipped case, flip its outputs back, and average. Returns the fused output matrix (`n_outputs × n_cases`), with the Mahalanobis penalty already applied to the confidence logit (row 1). """ -function neuralfoil_fused_output(params::KulfanParameters, alpha, Re; +function neuralfoil_fused_output(params, alpha, Re; model_size::String="xlarge", weights_dir=nothing, n_crit=9.0, xtr_upper=1.0, xtr_lower=1.0) - model = load_neuralfoil_model(model_size; weights_dir) x = prepare_inputs(params, alpha, Re; n_crit, xtr_upper, xtr_lower) + return fused_output(x, load_neuralfoil_model(model_size; weights_dir)) +end + +""" + fused_output(x, model) -> Matrix + +Symmetry-fused forward pass over a prepared input matrix, see +[`neuralfoil_fused_output`](@ref). Takes the inputs already built so a caller that +assembles its own batch does not go back through [`prepare_inputs`](@ref). +""" +function fused_output(x::AbstractMatrix, model::NeuralFoilModel) y = nn_forward(x, model) y[1, :] .-= squared_mahalanobis_distance(x, model) ./ (2 * model.n_inputs) diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index a1f3cdb9..f8ee8f7a 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -224,6 +224,13 @@ function calculate_stall_angle_list!(stall_angles::AbstractVector, # Default stall angle if none found panel_stall = stall_angle_if_none_detected + # A local expansion says nothing outside its fit window, so its curvature must + # not be read as a stall peak far from the operating point. + if panel.aero_model == TAYLOR + stall_angles[idx] = panel_stall + continue + end + # Start with minimum cl cl_old = cl_initial diff --git a/src/panel.jl b/src/panel.jl index 83edc848..aef44021 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -36,6 +36,7 @@ Represents a panel in a vortex step method simulation. All points and vectors ar ): Panel filaments, see: [BoundFilament](@ref) - `delta`::T=0: flap trailing-edge deflection [rad] - `crease_frac`::T=0: chordwise flap-hinge fraction (0–1); 0 disables the plate kink +- `alpha_ref`::Float64=0: expansion point [rad] of the `TAYLOR` coefficients """ @with_kw mutable struct Panel{T, CL, CD, CM, SA} TE_point_1::MVector{3, T} = zeros(MVector{3, T}) @@ -70,6 +71,7 @@ Represents a panel in a vortex step method simulation. All points and vectors ar ) delta::T = zero(T) crease_frac::T = zero(T) + alpha_ref::Float64 = 0.0 end """ @@ -216,6 +218,16 @@ function init_aero!(panel::Panel, section_1::Section, section_2::Section; panel.cl_coeffs = (c1[1] .+ c2[1]) ./ 2 panel.cd_coeffs = (c1[2] .+ c2[2]) ./ 2 panel.cm_coeffs = (c1[3] .+ c2[3]) ./ 2 + elseif panel.aero_model == TAYLOR + c1, c2 = section_1.aero_data, section_2.aero_data + (c1 isa Tuple{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}} && + c2 isa Tuple{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}}) || + throw(ArgumentError("TAYLOR requires aero_data = " * + "(alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs).")) + all(length.(c1[2:4]) .== length.(c2[2:4])) || + throw(ArgumentError("TAYLOR coefficient vectors must have equal length.")) + set_taylor_polar!(panel, (c1[1] + c2[1]) / 2, (c1[2] .+ c2[2]) ./ 2, + (c1[3] .+ c2[3]) ./ 2, (c1[4] .+ c2[4]) ./ 2) elseif !(panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES, INVISCID)) throw(ArgumentError("Unsupported aero model: $(panel.aero_model)")) end @@ -224,6 +236,29 @@ function init_aero!(panel::Panel, section_1::Section, section_2::Section; return nothing end +""" + set_taylor_polar!(panel, alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs) + +Overwrite a `TAYLOR` panel's local polar: the expansion point `alpha_ref` [rad] and the +ascending coefficients of the polynomials in `α - alpha_ref`. Copies into the panel's +existing coefficient vectors when the order is unchanged, so a live polar source can +refit every solve without allocating. The panel's aero model is set to `TAYLOR`. +""" +function set_taylor_polar!(panel::Panel, alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs) + panel.aero_model = TAYLOR + panel.alpha_ref = Float64(alpha_ref) + for (dst_sym, src) in ((:cl_coeffs, cl_coeffs), (:cd_coeffs, cd_coeffs), + (:cm_coeffs, cm_coeffs)) + dst = getfield(panel, dst_sym) + if length(dst) == length(src) + dst .= src + else + setfield!(panel, dst_sym, Vector{Float64}(src)) + end + end + return nothing +end + """ reinit!(panel, section_1, section_2, aero_center, control_point, bound_point_1, bound_point_2, x_airf, y_airf, z_airf, delta, vec; kwargs...) @@ -333,6 +368,8 @@ function calculate_cl(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} cl = 2 * cos(alpha) * sin(alpha)^2 end return R(cl) + elseif panel.aero_model == TAYLOR + return R(evalpoly(alpha - panel.alpha_ref, panel.cl_coeffs)) elseif panel.aero_model == INVISCID return R(2π * alpha) end @@ -361,6 +398,8 @@ function calculate_cd(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} return R(2 * sin(alpha)^3) end return R(evalpoly(rad2deg(alpha), panel.cd_coeffs)) + elseif panel.aero_model == TAYLOR + return R(evalpoly(alpha - panel.alpha_ref, panel.cd_coeffs)) elseif panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES) cd_interp = panel.cd_interp cd_interp === nothing && @@ -387,6 +426,8 @@ function calculate_cm(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} isnan(alpha) && return R(NaN) if panel.aero_model == POLY return R(evalpoly(rad2deg(alpha), panel.cm_coeffs)) + elseif panel.aero_model == TAYLOR + return R(evalpoly(alpha - panel.alpha_ref, panel.cm_coeffs)) elseif panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES) cm_interp = panel.cm_interp cm_interp === nothing && diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 1025ed41..3a4c5047 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -1287,6 +1287,20 @@ function calculate_new_aero_data(aero_model, return (alpha_left, delta_left, CL_data, CD_data, CM_data) + elseif isequal(model_type, TAYLOR) + data_left = aero_data[section_index] + data_right = aero_data[section_index + 1] + (data_left isa Tuple{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}} && + data_right isa Tuple{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}}) || + throw(ArgumentError("TAYLOR requires aero_data = " * + "(alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs).")) + return ( + data_left[1] * left_weight + data_right[1] * right_weight, + data_left[2] .* left_weight .+ data_right[2] .* right_weight, + data_left[3] .* left_weight .+ data_right[3] .* right_weight, + data_left[4] .* left_weight .+ data_right[4] .* right_weight, + ) + elseif isequal(model_type, POLY) data_left = aero_data[section_index] data_right = aero_data[section_index + 1] diff --git a/test/airfoil_aero/test_live_polar.jl b/test/airfoil_aero/test_live_polar.jl new file mode 100644 index 00000000..a867f5fc --- /dev/null +++ b/test/airfoil_aero/test_live_polar.jl @@ -0,0 +1,111 @@ +using Test +using LinearAlgebra +using Statistics +using VortexStepMethod +using VortexStepMethod.AirfoilAero +using VortexStepMethod: Panel, calculate_cl, calculate_cd, calculate_cm, + set_taylor_polar! + +@testset "TAYLOR panel polar" begin + panel = Panel{Float64}() + set_taylor_polar!(panel, deg2rad(5.0), [0.6, 5.0, -2.0], [0.02, 0.1, 1.0], + [-0.1, 0.2, 0.0]) + @test panel.aero_model == TAYLOR + alpha = deg2rad(6.0) + d = alpha - deg2rad(5.0) + @test calculate_cl(panel, alpha) ≈ 0.6 + 5.0d - 2.0d^2 + @test calculate_cd(panel, alpha) ≈ 0.02 + 0.1d + 1.0d^2 + @test calculate_cm(panel, alpha) ≈ -0.1 + 0.2d + + # A local expansion has no flap axis: delta must not change the answer. + @test calculate_cl(panel, alpha, deg2rad(9.0)) == calculate_cl(panel, alpha) + + coeffs = panel.cl_coeffs + set_taylor_polar!(panel, 0.0, [1.0, 2.0, 3.0], [1.0, 2.0, 3.0], [1.0, 2.0, 3.0]) + @test panel.cl_coeffs === coeffs # same order refits in place +end + +@testset "Kulfan deformation" begin + basis = KulfanBasis() + base = KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) + camber = @. 0.05 * basis.x^2 * (1 - basis.x) + deformed = deform_kulfan(basis, base, camber) + + shape = AirfoilAero.class_function(basis.x) .* + AirfoilAero.bernstein_basis(basis.x, 7) + for (weights, base_weights) in ((deformed.upper_weights, base.upper_weights), + (deformed.lower_weights, base.lower_weights)) + residual = shape * (weights .- base_weights) .- camber + @test maximum(abs, residual) < 0.01 * maximum(abs, camber) + end + @test deformed.leading_edge_weight == base.leading_edge_weight + @test deformed.TE_thickness == base.TE_thickness + + # Deforming the surfaces apart changes thickness, not just camber. + apart = deform_kulfan(basis, base, camber, -camber) + @test apart.upper_weights ≉ apart.lower_weights .+ + (base.upper_weights .- base.lower_weights) + + @test_throws ArgumentError deform_kulfan(basis, base, camber[1:end-1]) +end + +@testset "control point deflection" begin + basis = KulfanBasis() + sampled = control_point_deflection(basis, [0.0, 0.7, 0.8, 1.0], + [0.0, 0.03, 0.025, 0.0]) + @test sampled[1] ≈ 0.0 atol = 1e-12 + @test sampled[end] ≈ 0.0 atol = 1e-12 + @test maximum(sampled) ≈ 0.03 atol = 1e-3 + + # Unsorted input is the same curve as sorted input. + shuffled = control_point_deflection(basis, [0.8, 0.0, 1.0, 0.7], + [0.025, 0.0, 0.0, 0.03]) + @test shuffled ≈ sampled + @test_throws ArgumentError control_point_deflection(basis, [0.0, 0.0], [1.0, 2.0]) +end + +@testset "live polars" begin + base = KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) + n_panels = 3 + panels = [Panel{Float64}() for _ in 1:n_panels] + live = LivePolars(fill(base, n_panels)) + confidence = refresh_live_polars!(live, panels, deg2rad(6.0), 3e6) + @test 0.0 < confidence <= 1.0 + @test all(p -> p.aero_model == TAYLOR, panels) + + # The fit tracks a direct NeuralFoil sweep across its own window. + errors = [abs(calculate_cl(panels[1], deg2rad(6.0 + d)) - + neuralfoil_aero(base, 6.0 + d, 3e6).CL[1]) + for d in range(-4, 4, 17)] + @test maximum(errors) < 0.02 + + @test polar_drift(live, fill(deg2rad(6.0), n_panels)) ≈ 0.0 atol = 1e-12 + @test polar_drift(live, fill(deg2rad(10.0), n_panels)) ≈ 1.0 + + # Deforming the camber up raises lift at the same angle of attack. + flat = calculate_cl(panels[1], deg2rad(6.0)) + camber = @. 0.02 * live.basis.x * (1 - live.basis.x) + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6; + deflection=fill(camber, n_panels)) + @test calculate_cl(panels[1], deg2rad(6.0)) > flat + + @test_throws ArgumentError refresh_live_polars!(live, panels[1:2], 0.0, 3e6) + @test_throws ArgumentError LivePolars(fill(base, 2); + settings=LivePolarSettings(; order=4, n_samples=3)) +end + +@testset "TAYLOR spanwise blend" begin + left = (deg2rad(4.0), [0.5, 5.0, 0.0], [0.02, 0.0, 0.0], [-0.1, 0.0, 0.0]) + right = (deg2rad(8.0), [0.9, 5.0, 0.0], [0.06, 0.0, 0.0], [-0.3, 0.0, 0.0]) + blended = VortexStepMethod.calculate_new_aero_data( + (TAYLOR, TAYLOR), (left, right), 1, 0.25, 0.75) + @test blended[1] ≈ deg2rad(7.0) + @test blended[2] ≈ [0.8, 5.0, 0.0] + + section_left = Section([0.0, 1.0, 0.0], [1.0, 1.0, 0.0], TAYLOR, left) + section_right = Section([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], TAYLOR, right) + panel = Panel{Float64}() + VortexStepMethod.init_aero!(panel, section_left, section_right) + @test panel.alpha_ref ≈ deg2rad(6.0) + @test panel.cl_coeffs ≈ [0.7, 5.0, 0.0] +end diff --git a/test/runtests.jl b/test/runtests.jl index d613c2b5..08f8763c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -62,6 +62,7 @@ function include_selected_tests() should_run_test("wing_geometry/test_billowing.jl") && include("wing_geometry/test_billowing.jl") should_run_test("yaml_geometry/test_yaml_geometry.jl") && include("yaml_geometry/test_yaml_geometry.jl") should_run_test("airfoil_aero/test_airfoil_aero.jl") && include("airfoil_aero/test_airfoil_aero.jl") + should_run_test("airfoil_aero/test_live_polar.jl") && include("airfoil_aero/test_live_polar.jl") should_run_test("obj_adapter/test_obj_adapter.jl") && include("obj_adapter/test_obj_adapter.jl") should_run_test("surfplan/test_surfplan.jl") && include("surfplan/test_surfplan.jl") should_run_test("Aqua.jl") && include("Aqua.jl") From 26d90921c9feecf344172a9acf4acfc712538fc4 Mon Sep 17 00:00:00 2001 From: Bart Date: Tue, 25 Aug 2026 13:07:04 +0200 Subject: [PATCH 2/8] Continue a TAYLOR polar linearly past its fit window A quadratic fitted over +-4 deg is worse than useless ten degrees out, and that is exactly where a cold-started solve begins: with no previous angle of attack to expand about, the first fit sits at zero while the wing flies at fifteen. The arms then run away and the solve does not converge. Past the window, carry the polynomial's own value and slope at the edge on linearly instead. The fit is untouched inside its range, C1 at the edge, and the answer outside is an honest extrapolation the caller's drift guard then corrects by refitting where the solve went. Shrink-wrap a panel's contour before fitting its base Kulfan parameters, the route the offline polar generator already takes. Fitting a raw slice directly returns weights that oscillate by an order of magnitude more than any deformation will. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 5 ++- docs/src/private_functions.md | 1 + src/airfoil_aero/live_polar.jl | 19 +++++++---- src/panel.jl | 48 ++++++++++++++++++++++------ test/airfoil_aero/test_live_polar.jl | 21 ++++++++++++ 5 files changed, 78 insertions(+), 16 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9dac1907..147be68e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,7 +17,10 @@ NeuralFoil at a few angles per panel in one batched forward pass and least-squares fit each panel's `TAYLOR` polar, so a chordwise deformation reaches the aerodynamics as a shape change rather than a flap angle. - `polar_drift` reports how far the solve has left the fit window. + `polar_drift` reports how far the solve has left the fit window, and past that + window a `TAYLOR` polar is continued linearly off its edge value and slope + rather than following the polynomial's diverging arms, so a solve that steps + outside still converges and is corrected by the next refit. - `prepare_inputs` accepts one Kulfan shape per case, so a whole wing goes through NeuralFoil in a single forward pass. - Shared panel aerodynamics (`src/panel_aerodynamics.jl`): the per-panel physics diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 3453a80e..bb63e209 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -99,6 +99,7 @@ update_non_deformed_sections! ```@docs calculate_new_aero_data set_taylor_polar! +taylor_value assemble_polar_matrix load_matrix_polar_data read_aero_matrix diff --git a/src/airfoil_aero/live_polar.jl b/src/airfoil_aero/live_polar.jl index d1411c8c..f3a1ba2d 100644 --- a/src/airfoil_aero/live_polar.jl +++ b/src/airfoil_aero/live_polar.jl @@ -81,15 +81,21 @@ end panel_kulfan_parameters(panels; delta=0.0) -> Vector{KulfanParameters} Fit the undeformed Kulfan parameters of every panel from the surface contour its -`section_aero` carries, the starting point [`LivePolars`](@ref) deforms. Fitted once at -build time, never inside the loop. +`section_aero` carries, the starting point [`LivePolars`](@ref) deforms. + +The contour is [`shrink_wrap`](@ref)ped first, the same route the offline polar +generator takes: fitting a raw slice directly returns weights that oscillate by an +order of magnitude more than the deformation ever will, and neighbouring panels then +disagree enough to cost the VSM solve its convergence. The wrap costs ~15 ms a panel, +which is why this is done once at build time and never inside the loop. """ function panel_kulfan_parameters(panels; delta=0.0) return map(panels) do panel panel.section_aero === nothing && throw(ArgumentError( "Live polars need a section_aero contour on every panel.")) x, y, _, _ = section_surface(panel.section_aero, 0.0, delta) - fit_kulfan_parameters(collect(float.(x)), collect(float.(y))) + xw, yw = shrink_wrap(collect(float.(x)), collect(float.(y)), ShrinkWrap()) + return fit_kulfan_parameters(xw, yw, LeastSquaresFit()) end end @@ -150,9 +156,10 @@ function refresh_live_polars!(live::LivePolars, panels, alpha_ref, reynolds; cl, cd, cm, confidence = decode_coefficients(fused_output(live.inputs, model)) for i in 1:n_panels - window = ((i - 1) * n_samples + 1):(i * n_samples) - set_taylor_polar!(panels[i], alpha_vec[i], live.fit * cl[window], - live.fit * cd[window], live.fit * cm[window]) + samples = ((i - 1) * n_samples + 1):(i * n_samples) + set_taylor_polar!(panels[i], alpha_vec[i], live.fit * cl[samples], + live.fit * cd[samples], live.fit * cm[samples]; + window=live.settings.half_window) end return minimum(confidence) end diff --git a/src/panel.jl b/src/panel.jl index aef44021..13513d53 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -37,6 +37,7 @@ Represents a panel in a vortex step method simulation. All points and vectors ar - `delta`::T=0: flap trailing-edge deflection [rad] - `crease_frac`::T=0: chordwise flap-hinge fraction (0–1); 0 disables the plate kink - `alpha_ref`::Float64=0: expansion point [rad] of the `TAYLOR` coefficients +- `alpha_window`::Float64=0: half width [rad] the `TAYLOR` fit is valid over; 0 = unbounded """ @with_kw mutable struct Panel{T, CL, CD, CM, SA} TE_point_1::MVector{3, T} = zeros(MVector{3, T}) @@ -72,6 +73,7 @@ Represents a panel in a vortex step method simulation. All points and vectors ar delta::T = zero(T) crease_frac::T = zero(T) alpha_ref::Float64 = 0.0 + alpha_window::Float64 = 0.0 end """ @@ -237,16 +239,41 @@ function init_aero!(panel::Panel, section_1::Section, section_2::Section; end """ - set_taylor_polar!(panel, alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs) + taylor_value(coeffs, delta_alpha, window) -Overwrite a `TAYLOR` panel's local polar: the expansion point `alpha_ref` [rad] and the -ascending coefficients of the polynomials in `α - alpha_ref`. Copies into the panel's -existing coefficient vectors when the order is unchanged, so a live polar source can -refit every solve without allocating. The panel's aero model is set to `TAYLOR`. +A `TAYLOR` polar's coefficient at `delta_alpha = α - α_ref` [rad], continued linearly +beyond `±window`. A local expansion's arms diverge fast outside the range it was fitted +over — an order-2 fit is worse than useless a few degrees out — so past the edge the +polynomial's own value and slope there carry it on instead. That keeps a solve stepping +outside the window convergent and honest about being extrapolated, rather than chasing a +parabola to infinity; the caller's drift guard is what puts the fit back where the solve +went. `window = 0` leaves the polynomial unbounded. """ -function set_taylor_polar!(panel::Panel, alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs) +@inline function taylor_value(coeffs, delta_alpha, window) + if window > 0 && abs(delta_alpha) > window + edge = delta_alpha >= 0 ? window : -window + slope = sum((k - 1) * coeffs[k] * edge^(k - 2) for k in 2:length(coeffs); + init = 0.0) + return evalpoly(edge, coeffs) + slope * (delta_alpha - edge) + end + return evalpoly(delta_alpha, coeffs) +end + +""" + set_taylor_polar!(panel, alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs; window=0.0) + +Overwrite a `TAYLOR` panel's local polar: the expansion point `alpha_ref` [rad], the +ascending coefficients of the polynomials in `α - alpha_ref`, and the half width +`window` [rad] the fit is valid over, past which it is continued linearly (see +[`taylor_value`](@ref)). Copies into the panel's existing coefficient vectors when the +order is unchanged, so a live polar source can refit every solve without allocating. +The panel's aero model is set to `TAYLOR`. +""" +function set_taylor_polar!(panel::Panel, alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs; + window=0.0) panel.aero_model = TAYLOR panel.alpha_ref = Float64(alpha_ref) + panel.alpha_window = Float64(window) for (dst_sym, src) in ((:cl_coeffs, cl_coeffs), (:cd_coeffs, cd_coeffs), (:cm_coeffs, cm_coeffs)) dst = getfield(panel, dst_sym) @@ -369,7 +396,8 @@ function calculate_cl(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} end return R(cl) elseif panel.aero_model == TAYLOR - return R(evalpoly(alpha - panel.alpha_ref, panel.cl_coeffs)) + return R(taylor_value(panel.cl_coeffs, alpha - panel.alpha_ref, + panel.alpha_window)) elseif panel.aero_model == INVISCID return R(2π * alpha) end @@ -399,7 +427,8 @@ function calculate_cd(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} end return R(evalpoly(rad2deg(alpha), panel.cd_coeffs)) elseif panel.aero_model == TAYLOR - return R(evalpoly(alpha - panel.alpha_ref, panel.cd_coeffs)) + return R(taylor_value(panel.cd_coeffs, alpha - panel.alpha_ref, + panel.alpha_window)) elseif panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES) cd_interp = panel.cd_interp cd_interp === nothing && @@ -427,7 +456,8 @@ function calculate_cm(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} if panel.aero_model == POLY return R(evalpoly(rad2deg(alpha), panel.cm_coeffs)) elseif panel.aero_model == TAYLOR - return R(evalpoly(alpha - panel.alpha_ref, panel.cm_coeffs)) + return R(taylor_value(panel.cm_coeffs, alpha - panel.alpha_ref, + panel.alpha_window)) elseif panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES) cm_interp = panel.cm_interp cm_interp === nothing && diff --git a/test/airfoil_aero/test_live_polar.jl b/test/airfoil_aero/test_live_polar.jl index a867f5fc..0591cafb 100644 --- a/test/airfoil_aero/test_live_polar.jl +++ b/test/airfoil_aero/test_live_polar.jl @@ -25,6 +25,27 @@ using VortexStepMethod: Panel, calculate_cl, calculate_cd, calculate_cm, @test panel.cl_coeffs === coeffs # same order refits in place end +@testset "TAYLOR window is continued linearly" begin + window = deg2rad(4.0) + panel = Panel{Float64}() + set_taylor_polar!(panel, 0.0, [0.6, 5.0, -20.0], [0.02, 0.0, 0.0], + [-0.1, 0.0, 0.0]; window) + # Inside the window the polynomial is untouched. + @test calculate_cl(panel, 0.5window) ≈ 0.6 + 5.0 * 0.5window - 20.0 * (0.5window)^2 + # At the edge value and slope match, so the continuation is smooth. + edge = calculate_cl(panel, window) + slope = (calculate_cl(panel, window + 1e-7) - edge) / 1e-7 + inner = (edge - calculate_cl(panel, window - 1e-7)) / 1e-7 + @test slope ≈ inner rtol = 1e-4 + # And it stays linear rather than falling off with the parabola's arm, on both + # sides: the unbounded quadratic would be 0.6 + 5·d − 20·d² far out. + @test calculate_cl(panel, 4window) ≈ edge + slope * 3window rtol = 1e-5 + lower_edge = calculate_cl(panel, -window) + lower_slope = 5.0 - 2 * 20.0 * (-window) + @test calculate_cl(panel, -4window) ≈ lower_edge - lower_slope * 3window rtol = 1e-5 + @test calculate_cl(panel, 4window) > 0.6 + 5.0 * 4window - 20.0 * (4window)^2 +end + @testset "Kulfan deformation" begin basis = KulfanBasis() base = KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) From ec5c50e5f92cb231983415f74f7b55dc43959edf Mon Sep 17 00:00:00 2001 From: Bart Date: Tue, 25 Aug 2026 13:16:10 +0200 Subject: [PATCH 3/8] Say what a live polar refresh actually allocates It is not "nothing beyond the network's own forward pass": the deformation allocates a weight vector a panel and the deflection path allocates its resampled curve. Measured on a 44-panel wing, the deformation is 0.57 us and 608 B a panel and the forward pass is 3.6 ms and 6.2 MB, so the claim was wrong about the small part and right about which part matters. Co-Authored-By: Claude Opus 5 (1M context) --- src/airfoil_aero/live_polar.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/airfoil_aero/live_polar.jl b/src/airfoil_aero/live_polar.jl index f3a1ba2d..14b0a121 100644 --- a/src/airfoil_aero/live_polar.jl +++ b/src/airfoil_aero/live_polar.jl @@ -37,9 +37,10 @@ end LivePolars(base; settings=LivePolarSettings(), n_stations=60) Live polar source for a wing whose panels each carry one undeformed airfoil in `base`. -Holds the fixed CST basis, the per-panel expansion points and the scratch the refresh -reuses, so a steady-state refresh allocates nothing beyond the network's own forward -pass. Drive it with [`refresh_live_polars!`](@ref). +Holds the fixed CST basis, the per-panel expansion points and the network input scratch. +A refresh is dominated by the forward pass itself; the deformation is a matvec against +the constant basis, about half a microsecond a panel. Drive it with +[`refresh_live_polars!`](@ref). """ mutable struct LivePolars "Sampling and fitting configuration." From 7c325ec3b09f5c2d2f760b3851c748d468746d48 Mon Sep 17 00:00:00 2001 From: Bart Date: Tue, 25 Aug 2026 13:34:02 +0200 Subject: [PATCH 4/8] Draw the airfoil the live polars were actually generated from The lofted skin drew section_surface at the panel's delta, which under live polars is the undeformed tabulated contour: the picture showed the shape the solve was not flying. Deformation bugs could hide behind it. A panel now carries the deformed KulfanParameters its polar was generated from, written by the same call that installs the polar so the two cannot drift apart, and the skin lofts that. It is the object the airfoil solver was handed, not a second derivation of the same deformation, which is the point: if the deformation is wrong the drawing is wrong with it. KulfanParameters moves from AirfoilAero into the core module so a Panel can hold one; the fitting and coordinate routines stay where they were and both modules still export the type. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 6 +++++ docs/src/private_functions.md | 1 + docs/src/private_types.md | 2 +- ext/VortexStepMethodMakieExt.jl | 36 +++++++++++++++++++++------- src/VortexStepMethod.jl | 1 + src/airfoil_aero/AirfoilAero.jl | 3 ++- src/airfoil_aero/kulfan.jl | 18 -------------- src/airfoil_aero/live_polar.jl | 5 +++- src/panel.jl | 14 +++++++++-- src/section_aero.jl | 24 +++++++++++++++++++ test/airfoil_aero/test_live_polar.jl | 9 +++++++ 11 files changed, 87 insertions(+), 32 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 147be68e..af57bf31 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -23,6 +23,12 @@ outside still converges and is corrected by the next refit. - `prepare_inputs` accepts one Kulfan shape per case, so a whole wing goes through NeuralFoil in a single forward pass. +- A panel carries the deformed airfoil its polar was generated from as + `live_shape`, and the Makie lofted-airfoil skin draws that object rather than + the tabulated `delta` contour when it is there. It is the shape the solver was + handed, not a re-derivation, so a deformation bug is visible in the picture. + `KulfanParameters` moved from `AirfoilAero` into the core module for it; it is + still exported from both. - Shared panel aerodynamics (`src/panel_aerodynamics.jl`): the per-panel physics written once as pure, branch-free, number-type-generic functions of the section geometry and the flow — `panel_axes`, `panel_inflow`, `panel_force_directions`, diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index bb63e209..0691765a 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -226,6 +226,7 @@ map_airfoil_3d fitted_airfoil_3d generated_slices airfoil_skin_geometry +panel_contour panel_normal plate_hinge_local panel_plate_geometry diff --git a/docs/src/private_types.md b/docs/src/private_types.md index 21e223d3..ff63c1fc 100644 --- a/docs/src/private_types.md +++ b/docs/src/private_types.md @@ -7,6 +7,7 @@ CurrentModule = VortexStepMethod ### Wing Geometry, Panel and Aerodynamics ```@docs Panel +KulfanParameters PanelProperties Filament BoundFilament @@ -23,7 +24,6 @@ LEI_AIRFOIL_BREUKELS CurrentModule = VortexStepMethod.AirfoilAero ``` ```@docs -KulfanParameters KulfanBasis LivePolarSettings LivePolars diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl index 1356df73..7ec850d4 100644 --- a/ext/VortexStepMethodMakieExt.jl +++ b/ext/VortexStepMethodMakieExt.jl @@ -149,15 +149,34 @@ end """ airfoil_skin_geometry(body; R_b_w=nothing, T_b_w=nothing) -> (vertices, faces, ribs) -Lofted airfoil skin of a `BodyAerodynamics`: each section's deflected contour -(`section_surface` at the panel's `delta`) is fitted between the panel's `corner_points` -by a 2D similarity, TE pinned so a deflection bulges the fore body up. The skin reflects -`delta` only when the geometry carries per-`delta` slices (`obj_to_yaml` with a -`delta_range`); with δ=0-only data it renders undeflected, a deliberate cue that the -deflected slices are missing. Transformed to world by `R_b_w`/`T_b_w`. +Lofted airfoil skin of a `BodyAerodynamics`: each section's contour is fitted between +the panel's `corner_points` by a 2D similarity, TE pinned so a deflection bulges the +fore body up. + +The contour is the panel's own `live_shape` when it has one — the very +[`KulfanParameters`](@ref) object the live polar source deformed and handed to the +airfoil solver, not a re-derivation of it, so a deformation bug shows up in the picture +instead of being papered over. Otherwise it is `section_surface` at the panel's `delta`, +which reflects `delta` only when the geometry carries per-`delta` slices (`obj_to_yaml` +with a `delta_range`); with δ=0-only data it renders undeflected, a deliberate cue that +the deflected slices are missing. Transformed to world by `R_b_w`/`T_b_w`. `vertices`/`faces` triangulate the skin between consecutive equal-node sections; `ribs` is one closed contour polyline per section. Sections without contour data are skipped. """ +""" + panel_contour(panel, section) -> (x, y) + +The airfoil contour to draw a panel with: coordinates of its `live_shape` when a live +polar source has put one there, else the section's tabulated contour at the panel's +`delta`. The live branch reads the stored shape object itself, so what is drawn is what +was flown. +""" +function panel_contour(panel, section) + shape = panel.live_shape + isnothing(shape) || return VortexStepMethod.AirfoilAero.kulfan_to_coordinates(shape) + return VortexStepMethod.section_surface(section.section_aero, 0.0, panel.delta)[1:2] +end + function airfoil_skin_geometry(body; R_b_w=nothing, T_b_w=nothing) to_world(p) = (isnothing(R_b_w) || isnothing(T_b_w)) ? Point3f(p) : Point3f(R_b_w * p + T_b_w) @@ -169,10 +188,10 @@ function airfoil_skin_geometry(body; R_b_w=nothing, T_b_w=nothing) n_panels = n - 1 n_panels < 1 && continue for (i, section) in enumerate(sections) - isnothing(section.section_aero) && continue panel_idx = panel_offset + min(i, n_panels) panel_idx <= length(body.panels) || continue panel = body.panels[panel_idx] + (isnothing(section.section_aero) && isnothing(panel.live_shape)) && continue corners = panel.corner_points corner1 = Point3f(corners[:, 1]); corner3 = Point3f(corners[:, 3]) corner2 = Point3f(corners[:, 2]); corner4 = Point3f(corners[:, 4]) @@ -183,8 +202,7 @@ function airfoil_skin_geometry(body; R_b_w=nothing, T_b_w=nothing) chord_len = norm(chord) chord_len < 1e-9 && continue up = panel_normal(panel) - xs, ys, _, _ = VortexStepMethod.section_surface(section.section_aero, - 0.0, panel.delta) + xs, ys = panel_contour(panel, section) le_i = argmin(xs) le_x = xs[le_i]; le_y = ys[le_i] te_x = 0.5 * (xs[1] + xs[end]); te_y = 0.5 * (ys[1] + ys[end]) diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 93661842..57b0ebb9 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -37,6 +37,7 @@ export MVec3 export LLT, Model, VSM export AeroModel, INVISCID, POLY, LEI_AIRFOIL_BREUKELS, POLAR_MATRICES, POLAR_VECTORS, TAYLOR +export KulfanParameters export BILLOWING, COSINE, LINEAR, PanelDistribution, SPLIT_PROVIDED, UNCHANGED export ELLIPTIC, InitialGammaDistribution, ZEROS export FAILURE, FEASIBLE, INFEASIBLE, SolverStatus diff --git a/src/airfoil_aero/AirfoilAero.jl b/src/airfoil_aero/AirfoilAero.jl index 452f775a..524445e8 100644 --- a/src/airfoil_aero/AirfoilAero.jl +++ b/src/airfoil_aero/AirfoilAero.jl @@ -8,7 +8,8 @@ using NPZ using Xfoil using Printf: @sprintf using ..VortexStepMethod: SectionAero, interpolate_matrix_nans!, delta_suffix, - write_node_rows, section_surface, set_taylor_polar! + write_node_rows, section_surface, set_taylor_polar!, + KulfanParameters include("kulfan.jl") include("deform.jl") diff --git a/src/airfoil_aero/kulfan.jl b/src/airfoil_aero/kulfan.jl index aef6656d..dce90044 100644 --- a/src/airfoil_aero/kulfan.jl +++ b/src/airfoil_aero/kulfan.jl @@ -6,24 +6,6 @@ This module provides functions to: - Generate airfoil coordinates from Kulfan parameters """ -""" - KulfanParameters - -Kulfan CST parameters for an airfoil. - -# Fields -- `upper_weights::Vector{Float64}`: 8 weights for upper surface -- `lower_weights::Vector{Float64}`: 8 weights for lower surface -- `leading_edge_weight::Float64`: Leading edge modification weight -- `TE_thickness::Float64`: Trailing edge thickness -""" -struct KulfanParameters - upper_weights::Vector{Float64} - lower_weights::Vector{Float64} - leading_edge_weight::Float64 - TE_thickness::Float64 -end - """ KulfanFitMethod diff --git a/src/airfoil_aero/live_polar.jl b/src/airfoil_aero/live_polar.jl index 14b0a121..7edd7bae 100644 --- a/src/airfoil_aero/live_polar.jl +++ b/src/airfoil_aero/live_polar.jl @@ -124,6 +124,9 @@ deform the base airfoil by `deflection` (a chord-normalized deflection on `live.basis.x`, or `nothing` to keep the base shape), evaluate NeuralFoil at `alpha_ref ± half_window`, and least-squares fit the polynomial. +Each panel keeps the deformed shape it was evaluated at as its `live_shape`, so a plot +draws the airfoil the network actually saw. + `alpha_ref` [rad] and `reynolds` are per panel or scalar. Every panel and sample goes through the network in one forward pass. Returns the lowest analysis confidence over the batch, a value near zero meaning the deformed shape has left the region the network @@ -160,7 +163,7 @@ function refresh_live_polars!(live::LivePolars, panels, alpha_ref, reynolds; samples = ((i - 1) * n_samples + 1):(i * n_samples) set_taylor_polar!(panels[i], alpha_vec[i], live.fit * cl[samples], live.fit * cd[samples], live.fit * cm[samples]; - window=live.settings.half_window) + window=live.settings.half_window, shape=live.deformed[i]) end return minimum(confidence) end diff --git a/src/panel.jl b/src/panel.jl index 13513d53..34e830cc 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -38,6 +38,7 @@ Represents a panel in a vortex step method simulation. All points and vectors ar - `crease_frac`::T=0: chordwise flap-hinge fraction (0–1); 0 disables the plate kink - `alpha_ref`::Float64=0: expansion point [rad] of the `TAYLOR` coefficients - `alpha_window`::Float64=0: half width [rad] the `TAYLOR` fit is valid over; 0 = unbounded +- `live_shape`::Union{Nothing, KulfanParameters}=nothing: the deformed airfoil the polar was generated from """ @with_kw mutable struct Panel{T, CL, CD, CM, SA} TE_point_1::MVector{3, T} = zeros(MVector{3, T}) @@ -74,6 +75,7 @@ Represents a panel in a vortex step method simulation. All points and vectors ar crease_frac::T = zero(T) alpha_ref::Float64 = 0.0 alpha_window::Float64 = 0.0 + live_shape::Union{Nothing, KulfanParameters} = nothing end """ @@ -260,7 +262,8 @@ went. `window = 0` leaves the polynomial unbounded. end """ - set_taylor_polar!(panel, alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs; window=0.0) + set_taylor_polar!(panel, alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs; window=0.0, + shape=nothing) Overwrite a `TAYLOR` panel's local polar: the expansion point `alpha_ref` [rad], the ascending coefficients of the polynomials in `α - alpha_ref`, and the half width @@ -268,12 +271,19 @@ ascending coefficients of the polynomials in `α - alpha_ref`, and the half widt [`taylor_value`](@ref)). Copies into the panel's existing coefficient vectors when the order is unchanged, so a live polar source can refit every solve without allocating. The panel's aero model is set to `TAYLOR`. + +`shape` is the [`KulfanParameters`](@ref) the coefficients were generated from, stored +on the panel as `live_shape`. It is the object that was sampled, not a copy or a +re-derivation, and it is written here so a panel's polar and the shape behind it are +set together and cannot drift apart — which is what makes a plot of `live_shape` a +picture of what the solve actually flew rather than of what it should have. """ function set_taylor_polar!(panel::Panel, alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs; - window=0.0) + window=0.0, shape=nothing) panel.aero_model = TAYLOR panel.alpha_ref = Float64(alpha_ref) panel.alpha_window = Float64(window) + panel.live_shape = shape for (dst_sym, src) in ((:cl_coeffs, cl_coeffs), (:cd_coeffs, cd_coeffs), (:cm_coeffs, cm_coeffs)) dst = getfield(panel, dst_sym) diff --git a/src/section_aero.jl b/src/section_aero.jl index 3e219101..30dee3fd 100644 --- a/src/section_aero.jl +++ b/src/section_aero.jl @@ -1,3 +1,27 @@ +""" + KulfanParameters + +Kulfan CST parameters for an airfoil: the weights of the class-shape transformation +each surface is built from, the shared leading-edge modification weight and the +trailing-edge thickness. Lives here rather than in `AirfoilAero` because a +[`Panel`](@ref) carries the shape it is currently flying (`live_shape`), which a +plot or a traction pattern reads without knowing how the shape was produced. Fit one +with `AirfoilAero.fit_kulfan_parameters`, deform one with `AirfoilAero.deform_kulfan` +and turn one into coordinates with `AirfoilAero.kulfan_to_coordinates`. + +# Fields +- `upper_weights::Vector{Float64}`: weights for upper surface +- `lower_weights::Vector{Float64}`: weights for lower surface +- `leading_edge_weight::Float64`: Leading edge modification weight +- `TE_thickness::Float64`: Trailing edge thickness +""" +struct KulfanParameters + upper_weights::Vector{Float64} + lower_weights::Vector{Float64} + leading_edge_weight::Float64 + TE_thickness::Float64 +end + """ SectionAero diff --git a/test/airfoil_aero/test_live_polar.jl b/test/airfoil_aero/test_live_polar.jl index 0591cafb..785ca78d 100644 --- a/test/airfoil_aero/test_live_polar.jl +++ b/test/airfoil_aero/test_live_polar.jl @@ -110,6 +110,15 @@ end deflection=fill(camber, n_panels)) @test calculate_cl(panels[1], deg2rad(6.0)) > flat + # The panel keeps the object that was sampled, not a copy of it: a plot of + # live_shape has to be a picture of what the solve flew. + @test panels[1].live_shape === live.deformed[1] + @test panels[1].live_shape.upper_weights != + live.base[1].upper_weights # the deformation reached it + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6) + @test panels[1].live_shape === live.deformed[1] + @test panels[1].live_shape.upper_weights ≈ live.base[1].upper_weights + @test_throws ArgumentError refresh_live_polars!(live, panels[1:2], 0.0, 3e6) @test_throws ArgumentError LivePolars(fill(base, 2); settings=LivePolarSettings(; order=4, n_samples=3)) From d72601b2192dd368cb5acd06200f7a4b6a8fd80a Mon Sep 17 00:00:00 2001 From: Bart Date: Tue, 25 Aug 2026 13:39:47 +0200 Subject: [PATCH 5/8] Put panel_contour ahead of the skin's docstring The new docstring landed between airfoil_skin_geometry's docstring and its function, so the first was documenting a string literal and the whole Makie extension failed to precompile: not just airfoil skins, every plot VortexStepMethod draws. It only surfaced on loading the extension for real, which is the argument for the test this adds. The test puts a cambered live_shape on every panel and checks the rib node count switches from the tabulated contour's to the Kulfan one's and back when it is cleared, asserting first that the two counts differ so it cannot pass by drawing the wrong thing. Co-Authored-By: Claude Opus 5 (1M context) --- ext/VortexStepMethodMakieExt.jl | 28 ++++++++++++++-------------- test/plotting/test_plotting.jl | 22 ++++++++++++++++++++++ 2 files changed, 36 insertions(+), 14 deletions(-) diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl index 7ec850d4..ddd010e2 100644 --- a/ext/VortexStepMethodMakieExt.jl +++ b/ext/VortexStepMethodMakieExt.jl @@ -146,6 +146,20 @@ function Makie.plot!(ax, panel::VortexStepMethod.Panel; color=(:red, 0.2), R_b_w return plots end +""" + panel_contour(panel, section) -> (x, y) + +The airfoil contour to draw a panel with: coordinates of its `live_shape` when a live +polar source has put one there, else the section's tabulated contour at the panel's +`delta`. The live branch reads the stored shape object itself, so what is drawn is what +was flown. +""" +function panel_contour(panel, section) + shape = panel.live_shape + isnothing(shape) || return VortexStepMethod.AirfoilAero.kulfan_to_coordinates(shape) + return VortexStepMethod.section_surface(section.section_aero, 0.0, panel.delta)[1:2] +end + """ airfoil_skin_geometry(body; R_b_w=nothing, T_b_w=nothing) -> (vertices, faces, ribs) @@ -163,20 +177,6 @@ the deflected slices are missing. Transformed to world by `R_b_w`/`T_b_w`. `vertices`/`faces` triangulate the skin between consecutive equal-node sections; `ribs` is one closed contour polyline per section. Sections without contour data are skipped. """ -""" - panel_contour(panel, section) -> (x, y) - -The airfoil contour to draw a panel with: coordinates of its `live_shape` when a live -polar source has put one there, else the section's tabulated contour at the panel's -`delta`. The live branch reads the stored shape object itself, so what is drawn is what -was flown. -""" -function panel_contour(panel, section) - shape = panel.live_shape - isnothing(shape) || return VortexStepMethod.AirfoilAero.kulfan_to_coordinates(shape) - return VortexStepMethod.section_surface(section.section_aero, 0.0, panel.delta)[1:2] -end - function airfoil_skin_geometry(body; R_b_w=nothing, T_b_w=nothing) to_world(p) = (isnothing(R_b_w) || isnothing(T_b_w)) ? Point3f(p) : Point3f(R_b_w * p + T_b_w) diff --git a/test/plotting/test_plotting.jl b/test/plotting/test_plotting.jl index 14807a6f..e3c5dbc8 100644 --- a/test/plotting/test_plotting.jl +++ b/test/plotting/test_plotting.jl @@ -521,6 +521,28 @@ end # In-place update refreshes the registered skin observables without error. @test_nowarn Makie.plot!(body_aero) + # Live polars: the skin has to draw the panel's stored deformed shape, not the + # tabulated contour, or a deformation bug is invisible in the picture. + basis = VortexStepMethod.AirfoilAero.KulfanBasis() + base = VortexStepMethod.KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) + cambered = VortexStepMethod.AirfoilAero.deform_kulfan(basis, base, + @. 0.08 * basis.x * (1 - basis.x)) + for panel in body_aero.panels + VortexStepMethod.set_taylor_polar!(panel, 0.0, [0.5, 5.0, 0.0], + [0.02, 0.0, 0.0], [-0.05, 0.0, 0.0]; window=deg2rad(4), shape=cambered) + end + @test body_aero.panels[1].live_shape === cambered + v_live, f_live, ribs_live = airfoil_skin_geometry(body_aero) + live_nodes = length(VortexStepMethod.AirfoilAero.kulfan_to_coordinates(cambered)[1]) + @test live_nodes != n_node # the two routes are distinguishable + @test all(rib -> length(rib) == live_nodes, ribs_live) + @test length(f_live) == 2 * (n_sections - 1) * (live_nodes - 1) + @test_nowarn Makie.plot!(Axis3(Figure()[1, 1]), body_aero; airfoils=true) + for panel in body_aero.panels + panel.live_shape = nothing + end + @test all(rib -> length(rib) == n_node, airfoil_skin_geometry(body_aero)[3]) + # A body without any surface aero yields no skin geometry, but still plots. plain_body = create_body_aero() v_plain, f_plain, ribs_plain = airfoil_skin_geometry(plain_body) From 9ef29f87b785a2a012f3a3300b31873e04a5aec0 Mon Sep 17 00:00:00 2001 From: Bart Date: Tue, 25 Aug 2026 14:19:23 +0200 Subject: [PATCH 6/8] Take the chord line out of a deflection, and never return a negative drag Two things the SK100 found. C(x) vanishes at both chord ends, so the CST basis can put neither edge off the chord line: a deflection that does not end at zero is asking for a chord rotation and translation. Projecting it anyway is not merely inexact, it is unstable, and pinv answered a 2.5% chord deflection with weights three times larger than the camber they were meant to carry. chord_residual removes the line through the deflection's own endpoints; chord_line hands it back, since it is a chord rotation the caller may still owe its angle of attack. A drag fit has its minimum inside the window, so the slope at the lower edge points down and the linear continuation carried it through zero. A negative drag coefficient is an energy source, not an extrapolation; calculate_cd floors it. Neither of these was what broke the SK100 settle, which was a bring-up transient slewing alpha twenty degrees out of the window in ten milliseconds. Both are wrong on their own terms. Co-Authored-By: Claude Opus 5 (1M context) --- docs/src/functions.md | 2 ++ src/airfoil_aero/AirfoilAero.jl | 1 + src/airfoil_aero/deform.jl | 45 ++++++++++++++++++++++++++-- src/panel.jl | 12 ++++++-- test/airfoil_aero/test_live_polar.jl | 44 +++++++++++++++++++++++++++ 5 files changed, 99 insertions(+), 5 deletions(-) diff --git a/docs/src/functions.md b/docs/src/functions.md index 380b8523..94b45af4 100644 --- a/docs/src/functions.md +++ b/docs/src/functions.md @@ -42,6 +42,8 @@ analyze_section analyze_sweep neuralfoil_aero deform_kulfan +chord_residual +chord_line control_point_deflection panel_kulfan_parameters refresh_live_polars! diff --git a/src/airfoil_aero/AirfoilAero.jl b/src/airfoil_aero/AirfoilAero.jl index 524445e8..abf345dc 100644 --- a/src/airfoil_aero/AirfoilAero.jl +++ b/src/airfoil_aero/AirfoilAero.jl @@ -33,6 +33,7 @@ export fit_kulfan_parameters, kulfan_to_coordinates export NeuralFoilModel, NeuralFoilResult, load_neuralfoil_model export neuralfoil_aero, neuralfoil_section export KulfanBasis, deform_kulfan, control_point_deflection +export chord_residual, chord_line export LivePolarSettings, LivePolars, panel_kulfan_parameters export refresh_live_polars!, polar_drift export AbstractAirfoilSolver, XFoilSolver, NeuralFoilSolver diff --git a/src/airfoil_aero/deform.jl b/src/airfoil_aero/deform.jl index 52a165f5..08b61779 100644 --- a/src/airfoil_aero/deform.jl +++ b/src/airfoil_aero/deform.jl @@ -43,6 +43,12 @@ deflection to both surfaces, leaving the thickness distribution untouched; the four-argument form deforms the surfaces independently, which is what a double-skin membrane with its own upper and lower control points needs. +Each deflection is first reduced to the part the basis can carry, see +[`chord_residual`](@ref): the straight line through its own endpoints is a chord +rotation and translation, which the basis cannot express and `pinv` answers with +runaway weights. Take that line from [`chord_line`](@ref) if the frame the deflection +was measured in has not already absorbed it. + The leading-edge weight and the trailing-edge thickness are carried over unchanged: they are the two shape freedoms a chord-referenced deflection cannot resolve. """ @@ -55,9 +61,42 @@ function deform_kulfan(basis::KulfanBasis, base::KulfanParameters, (length(upper_deflection) == length(basis.x) && length(lower_deflection) == length(basis.x)) || throw(ArgumentError( "Deflections must be sampled on the basis' $(length(basis.x)) stations.")) - return KulfanParameters(base.upper_weights .+ basis.projection * upper_deflection, - base.lower_weights .+ basis.projection * lower_deflection, - base.leading_edge_weight, base.TE_thickness) + return KulfanParameters( + base.upper_weights .+ basis.projection * chord_residual(basis, upper_deflection), + base.lower_weights .+ basis.projection * chord_residual(basis, lower_deflection), + base.leading_edge_weight, base.TE_thickness) +end + +""" + chord_residual(basis, deflection) -> Vector{Float64} + chord_line(basis, deflection) -> (offset, slope) + +The part of a deflection the CST basis can represent, and the part it cannot. `C(x)` +vanishes at both chord ends, so the basis can put neither the leading nor the trailing +edge off the chord line: a deflection that does not end at zero is asking for a chord +**rotation and translation**, not a camber change. Projecting it anyway is not merely +inexact, it is unstable — `pinv` answers an unrepresentable end displacement with +weights an order of magnitude past the ones it is correcting, and the airfoil that +comes back is not one. + +[`chord_residual`](@ref) removes the straight line through the deflection's own +endpoints and returns what is left. [`chord_line`](@ref) returns that line as +`(offset, slope)`, both over chord: `offset` moves the leading edge and `slope` is the +chord rotation `atan(slope)` [rad] the caller owes its angle of attack, when the frame +it measured the deflection in has not already absorbed it. +""" +function chord_residual(basis::KulfanBasis, deflection::AbstractVector) + offset, slope = chord_line(basis, deflection) + return deflection .- (offset .+ slope .* basis.x) +end + +function chord_line(basis::KulfanBasis, deflection::AbstractVector) + length(deflection) == length(basis.x) || throw(ArgumentError( + "Deflection must be sampled on the basis' $(length(basis.x)) stations.")) + span = basis.x[end] - basis.x[1] + abs(span) < eps() && return (deflection[1], 0.0) + slope = (deflection[end] - deflection[1]) / span + return (deflection[1] - slope * basis.x[1], slope) end deform_kulfan(basis::KulfanBasis, base::KulfanParameters, diff --git a/src/panel.jl b/src/panel.jl index 34e830cc..e29d1f73 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -250,6 +250,11 @@ polynomial's own value and slope there carry it on instead. That keeps a solve s outside the window convergent and honest about being extrapolated, rather than chasing a parabola to infinity; the caller's drift guard is what puts the fit back where the solve went. `window = 0` leaves the polynomial unbounded. + +The continuation is a straight line, so it says nothing about stall and can carry a +coefficient somewhere physically impossible. [`calculate_cd`](@ref) therefore floors +its result at zero: every other coefficient may be extrapolated, but a negative drag +would feed energy into whatever reads it. """ @inline function taylor_value(coeffs, delta_alpha, window) if window > 0 && abs(delta_alpha) > window @@ -437,8 +442,11 @@ function calculate_cd(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} end return R(evalpoly(rad2deg(alpha), panel.cd_coeffs)) elseif panel.aero_model == TAYLOR - return R(taylor_value(panel.cd_coeffs, alpha - panel.alpha_ref, - panel.alpha_window)) + # A drag fit has its minimum inside the window, so the slope at the lower + # edge points down and the continuation would carry it through zero into a + # negative drag — an energy source, not an extrapolation. + return R(max(zero(R), taylor_value(panel.cd_coeffs, alpha - panel.alpha_ref, + panel.alpha_window))) elseif panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES) cd_interp = panel.cd_interp cd_interp === nothing && diff --git a/test/airfoil_aero/test_live_polar.jl b/test/airfoil_aero/test_live_polar.jl index 785ca78d..5b443bf6 100644 --- a/test/airfoil_aero/test_live_polar.jl +++ b/test/airfoil_aero/test_live_polar.jl @@ -37,6 +37,20 @@ end slope = (calculate_cl(panel, window + 1e-7) - edge) / 1e-7 inner = (edge - calculate_cl(panel, window - 1e-7)) / 1e-7 @test slope ≈ inner rtol = 1e-4 + # Drag is the exception: its fit has a minimum inside the window, so continuing + # the lower edge's slope would take it through zero. A negative drag coefficient + # is an energy source, and on the SK100 it drove the apparent wind from 12 to + # 21.8 m/s inside one step before the solve died. + drag = Panel{Float64}() + set_taylor_polar!(drag, 0.0, [0.5, 0.0, 0.0], [0.02, 0.4, 6.0], + [-0.05, 0.0, 0.0]; window) + @test calculate_cd(drag, -0.5window) < calculate_cd(drag, 0.0) # still falling + @test calculate_cd(drag, -window) >= 0.0 + @test calculate_cd(drag, -4window) >= 0.0 + @test calculate_cd(drag, -20window) >= 0.0 + # Above the window the continuation is untouched, drag grows. + @test calculate_cd(drag, 4window) > calculate_cd(drag, window) > 0.0 + # And it stays linear rather than falling off with the parabola's arm, on both # sides: the unbounded quadratic would be 0.6 + 5·d − 20·d² far out. @test calculate_cl(panel, 4window) ≈ edge + slope * 3window rtol = 1e-5 @@ -70,6 +84,36 @@ end @test_throws ArgumentError deform_kulfan(basis, base, camber[1:end-1]) end +@testset "an unrepresentable chord line is removed, not projected" begin + basis = KulfanBasis() + base = KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) + camber = @. 0.02 * basis.x * (1 - basis.x) + + # A deflection whose ends are off zero is a chord rotation and translation. The + # basis cannot hold either, and projecting one anyway used to answer a 2.5% chord + # deflection with weights four times the ones it was correcting. + tilted = camber .+ 0.03 .* basis.x .+ 0.01 + offset, slope = chord_line(basis, tilted) + @test offset ≈ 0.01 atol = 1e-12 + @test slope ≈ 0.03 atol = 1e-12 + residual = chord_residual(basis, tilted) + @test residual ≈ camber atol = 1e-12 + @test abs(residual[1]) < 1e-12 && abs(residual[end]) < 1e-12 + + # So the tilted deflection has to give the same airfoil as the pure camber, and + # weights of the order of the camber rather than of the tilt. + pure = deform_kulfan(basis, base, camber) + tilt = deform_kulfan(basis, base, tilted) + @test tilt.upper_weights ≈ pure.upper_weights + @test maximum(abs, tilt.upper_weights .- base.upper_weights) < + 10 * maximum(abs, camber) + + # A pure chord rotation is no shape change at all. + rotation = deform_kulfan(basis, base, 0.05 .* basis.x) + @test rotation.upper_weights ≈ base.upper_weights atol = 1e-12 + @test rotation.lower_weights ≈ base.lower_weights atol = 1e-12 +end + @testset "control point deflection" begin basis = KulfanBasis() sampled = control_point_deflection(basis, [0.0, 0.7, 0.8, 1.0], From d7838bce001832cb638166de1530e0216f656380 Mon Sep 17 00:00:00 2001 From: Bart Date: Tue, 25 Aug 2026 14:30:51 +0200 Subject: [PATCH 7/8] Keep generated aero tables out of git VortexStepMethod writes .arrow node tables itself, so generating polars inside the repo leaves them stageable. Nothing has been committed, and the two packages downstream already ignore them. Co-Authored-By: Claude Opus 5 (1M context) --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitignore b/.gitignore index d7c105ab..16f91689 100644 --- a/.gitignore +++ b/.gitignore @@ -59,3 +59,6 @@ lib/**/Manifest.toml # Generated test geometries (regenerated on demand by ram_air_matrix_wing) test/generated/ + +# Generated aero tables and run logs +*.arrow From 6c0ccbe381be6c1a55c6cdd760e293424d05812b Mon Sep 17 00:00:00 2001 From: Bart Date: Wed, 26 Aug 2026 15:03:27 +0200 Subject: [PATCH 8/8] Sample a live polar instead of fitting one MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A local polynomial cannot hold a stall. Measured on the SK100 over a causal sweep, no order/window combination the fit offered ever turned CL over: the true polar peaks at 24.0 deg with a lift slope collapsing to 0.0044/deg, every fit reported 8x to 13x that, and narrowing the window made it worse because the fit is centred on the previous step's alpha and never reaches the knee the solve is walking into. On a saved failing state that mechanism killed the run — one tip panel's window swallowed the post-peak collapse, its slope flipped to -0.0230/deg where the local slope was +0.0503/deg, and the solve went non-finite. SAMPLED replaces TAYLOR: cl/cd/cm at ascending angles per panel, linear between them and the end value held past either end. The samples are the polar, so a knee between two of them is represented rather than averaged into a slope, and past the ends the answer is stale rather than wrong-signed. The knots move with the panel and are written in place, so a refresh costs the network pass and nothing else. TAYLOR is removed rather than kept: it never shipped, and two local-polar models where one now works is one too many. Live surface pressure follows the same shapes. refresh_live_pressure! takes Cp from one batched pass at the converged alpha, contour_pressure resamples it onto the panel's own contour nodes, and live_surface_friction! gives skin friction from the flat-plate closure at the panel's live Reynolds. contour_shape_matrix and live_shape_offset! then move the contour itself by the deformation's own camber increment — a traction is -Cp*n with n a finite difference of neighbouring node positions, so an undeformed contour points the whole load the wrong way wherever the section has moved. The chord fractions never move, only the offsets, so a surface-to-point map built on that contour stays valid. decode_surface_pressure and contour_pressure also back the NeuralFoilSolver contour assembly, so the live path and the offline table generator are one implementation; a test pins them equal at rtol 1e-8. Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 400 +++++++++++------- docs/src/functions.md | 6 + docs/src/private_functions.md | 6 +- src/VortexStepMethod.jl | 26 +- src/airfoil_aero/AirfoilAero.jl | 6 +- .../airfoil_solvers/neuralfoil_solver.jl | 5 +- src/airfoil_aero/deform.jl | 26 +- src/airfoil_aero/live_polar.jl | 254 ++++++++--- src/airfoil_aero/neuralfoil.jl | 40 +- src/body_aerodynamics.jl | 6 +- src/panel.jl | 123 +++--- src/wing_geometry.jl | 14 - test/airfoil_aero/test_live_polar.jl | 228 ++++++---- test/plotting/test_plotting.jl | 5 +- 14 files changed, 736 insertions(+), 409 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index af57bf31..5a97726d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,24 +3,38 @@ ## VortexStepMethod v4.2.0 2026-08-25 ### Added -- `TAYLOR` aero model: a per-panel polynomial of arbitrary order in `α - α_ref` - [rad], the local expansion a live polar source refits every solve. Sections - carry `aero_data = (alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs)` and blend - spanwise like any other model; `set_taylor_polar!` rewrites a panel's fit in - place. It is valid only inside its fit window, so it ignores `delta` and is - skipped by the stall-angle scan. + +- `SAMPLED` aero model: cl/cd/cm sampled at ascending angles of attack per + panel, interpolated between them and held flat past either end. This is what a + live polar source writes every solve; `set_sampled_polar!` rewrites a panel's + knots and values in place. It spans only the angles it was sampled over, so it + ignores `delta` and is skipped by the stall-angle scan — but inside that range + it represents a stall, which is what a local fit cannot do. - Live in-memory polars in `AirfoilAero`: `KulfanBasis` and `deform_kulfan` deform a fixed Kulfan fit analytically (a matvec against a constant CST basis, never a refit, which is non-unique); `control_point_deflection` resamples a deflection given at arbitrary chord fractions — beam nodes, membrane nodes — onto that basis. `LivePolars` and `refresh_live_polars!` then evaluate - NeuralFoil at a few angles per panel in one batched forward pass and - least-squares fit each panel's `TAYLOR` polar, so a chordwise deformation - reaches the aerodynamics as a shape change rather than a flap angle. - `polar_drift` reports how far the solve has left the fit window, and past that - window a `TAYLOR` polar is continued linearly off its edge value and slope - rather than following the polynomial's diverging arms, so a solve that steps - outside still converges and is corrected by the next refit. + NeuralFoil on a grid of angles per panel in one batched forward pass and write + those values in as each panel's `SAMPLED` polar, so a chordwise deformation + reaches the aerodynamics as a shape change rather than a flap angle. The grid + moves with the panel — `LivePolarSettings` carries the offsets off its current + angle of attack — and `polar_drift` reports how far the solve has left it, + past which the polar holds its last sampled value instead of extrapolating. +- Live surface pressure in `AirfoilAero`: `refresh_live_pressure!` regenerates + every panel's `Cp` from its deformed shape in one batched forward pass at the + converged angle of attack, resampled onto the panel's own contour nodes by + `contour_pressure`; `live_surface_friction!` gives the matching skin friction + from the flat-plate closure at the panel's live Reynolds. A caller that + spreads panel forces over a surface can now make the pattern track the + deformation the forces already do. `contour_shape_matrix` and + `live_shape_offset!` move the contour itself by the same camber increment the + airfoil was deformed by, so the normals and segment areas a traction is built + from are the deformed section's rather than the reference one's — the chord + fractions never move, only the offsets, so a surface→point map built on that + contour stays valid. `decode_surface_pressure` and `contour_pressure` also + back the `NeuralFoilSolver` contour assembly, so the live path and the offline + table generator are one implementation. - `prepare_inputs` accepts one Kulfan shape per case, so a whole wing goes through NeuralFoil in a single forward pass. - A panel carries the deformed airfoil its polar was generated from as @@ -30,22 +44,24 @@ `KulfanParameters` moved from `AirfoilAero` into the core module for it; it is still exported from both. - Shared panel aerodynamics (`src/panel_aerodynamics.jl`): the per-panel physics - written once as pure, branch-free, number-type-generic functions of the section - geometry and the flow — `panel_axes`, `panel_inflow`, `panel_force_directions`, - `panel_loads` and the small helpers around them. `update_panel_properties!`, - `init_pos!`, `calc_forces!` and `calculate_results` now call them instead of - spelling the algebra out three times, and `SymbolicAWEModels` traces the same - functions with symbolic arguments to build its equations, so the two packages - can no longer drift apart. Panel geometry is unchanged bit for bit; forces, - moments and coefficients agree to within 1 ulp, the products having been - reassociated. `calc_forces!` stays zero-allocation. + written once as pure, branch-free, number-type-generic functions of the + section geometry and the flow — `panel_axes`, `panel_inflow`, + `panel_force_directions`, `panel_loads` and the small helpers around them. + `update_panel_properties!`, `init_pos!`, `calc_forces!` and + `calculate_results` now call them instead of spelling the algebra out three + times, and `SymbolicAWEModels` traces the same functions with symbolic + arguments to build its equations, so the two packages can no longer drift + apart. Panel geometry is unchanged bit for bit; forces, moments and + coefficients agree to within 1 ulp, the products having been reassociated. + `calc_forces!` stays zero-allocation. - `effective_alpha` and the `deficiency` argument of `panel_inflow` carry an unsteady lag (a Wagner indicial deficiency) into the angle the polars are read at, while the geometric angle still turns the force. Unused by the solver, - which has no unsteady state; it is shared so a symbolic consumer that does have - one reads the same definition. + which has no unsteady state; it is shared so a symbolic consumer that does + have one reads the same definition. ### Changed + - `section_pitch_rate` gained a three-argument form taking the trailing minus leading edge apparent wind directly. The four-argument form is unchanged. - Norms inside the shared panel aerodynamics are floored (`smooth_norm`) rather @@ -55,37 +71,43 @@ ## VortexStepMethod v4.1.2 2026-08-22 ### Fixed + - Billowing now takes its rotation axis from the section pair itself instead of projecting the leading-edge span vector on `spanwise_direction`. Sections run - `+y` to `-y`, so the pair already fixes the sign. The projection was decided by - 3 % of the vector on a C-shaped kite's tip panels — 11 mm of span between the - outermost sections — so a deforming tip could flip that one panel's billow. + `+y` to `-y`, so the pair already fixes the sign. The projection was decided + by 3 % of the vector on a C-shaped kite's tip panels — 11 mm of span between + the outermost sections — so a deforming tip could flip that one panel's + billow. ## VortexStepMethod v4.1.1 2026-08-19 ### Fixed + - Section twist now rotates the chord fully about the spanwise axis. The axial Rodrigues term was missing, so a swept or dihedral section's along-span chord - component was scaled by `cos(theta)`: the chord shortened and tilted out of the - twist plane, growing as `theta^2`. + component was scaled by `cos(theta)`: the chord shortened and tilted out of + the twist plane, growing as `theta^2`. ## VortexStepMethod v4.1.0 2026-08-17 ### Added -- `table_format` keyword on `write_section_aero`, `generate_airfoils`, `obj_to_yaml` - and `surfplan_to_aero_yaml`: `:csv` (default, readable) or `:arrow` (binary, ~40× - faster to load and 2.6× smaller). `read_section_aero` detects the format from the - file suffix, so a geometry YAML can reference either. + +- `table_format` keyword on `write_section_aero`, `generate_airfoils`, + `obj_to_yaml` and `surfplan_to_aero_yaml`: `:csv` (default, readable) or + `:arrow` (binary, ~40× faster to load and 2.6× smaller). `read_section_aero` + detects the format from the file suffix, so a geometry YAML can reference + either. - `geometry_path` keyword on `obj_to_yaml`, naming the geometry YAML itself instead of always writing `output_dir/geometry.yaml`. Point it outside the - table directory and the emitted table references carry the path from the YAML's - directory to `output_dir`, which is what the geometry loader resolves them - against — so a generated dataset can keep its bulk in a subdirectory while the - geometry sits with the hand-written ones. -- `convert_node_table` and `write_node_rows` rewrite a per-node table in the format - the destination suffix names. `obj_to_yaml` migrates an existing dataset with - them when `table_format` differs from what the directory holds, so a dataset - changes format without re-running the airfoil solver that produced it. + table directory and the emitted table references carry the path from the + YAML's directory to `output_dir`, which is what the geometry loader resolves + them against — so a generated dataset can keep its bulk in a subdirectory + while the geometry sits with the hand-written ones. +- `convert_node_table` and `write_node_rows` rewrite a per-node table in the + format the destination suffix names. `obj_to_yaml` migrates an existing + dataset with them when `table_format` differs from what the directory holds, + so a dataset changes format without re-running the airfoil solver that + produced it. - `flow_curvature` solver setting (default `false`). When enabled, each section gets the thin-airfoil pitch-rate moment increment `Δcm = -(π/4) q̂` with `q̂ = q c / (2 v_rel)`. A section rotating about its own spanwise axis sees an @@ -94,131 +116,146 @@ point cannot represent. The lift response to `q` was already exact because the inflow is sampled at the three-quarter-chord point, so only the moment was missing. -- `pitch_rate_dist` field on `BodyAerodynamics`, holding the `q` above per panel. - `set_va!(body_aero, va, omega)` fills it by projecting `omega` onto every - panel's `y_airf`, so panels at different dihedral see different rates from one - body rate. The distributed `set_va!(body_aero, va_distribution; - pitch_rate_dist)` takes it directly, so twist and flapping rates of a deforming - wing — which no single body rate can express — reach the moment. Omitting the - keyword zeroes it rather than reusing a stale `omega`, which the distributed - form never sets. -- `section_pitch_rate(velocity_leading, velocity_trailing, z_airf, chord)` builds - one entry of that distribution from a section's edge velocities, and reduces to - `ω ⋅ y_airf` for rigid motion. -- `bin/update_default_manifests` regenerates both checked-in default manifests in - one command, and `bin/install` takes `+X.Y` / `--version X.Y` to pick a Julia - channel without prompting. Selecting a version no longer moves the juliaup - default, so regenerating the 1.11 manifest leaves the shell on whatever channel - it was using. +- `pitch_rate_dist` field on `BodyAerodynamics`, holding the `q` above per + panel. `set_va!(body_aero, va, omega)` fills it by projecting `omega` onto + every panel's `y_airf`, so panels at different dihedral see different rates + from one body rate. The distributed + `set_va!(body_aero, va_distribution; + pitch_rate_dist)` takes it directly, so + twist and flapping rates of a deforming wing — which no single body rate can + express — reach the moment. Omitting the keyword zeroes it rather than reusing + a stale `omega`, which the distributed form never sets. +- `section_pitch_rate(velocity_leading, velocity_trailing, z_airf, chord)` + builds one entry of that distribution from a section's edge velocities, and + reduces to `ω ⋅ y_airf` for rigid motion. +- `bin/update_default_manifests` regenerates both checked-in default manifests + in one command, and `bin/install` takes `+X.Y` / `--version X.Y` to pick a + Julia channel without prompting. Selecting a version no longer moves the + juliaup default, so regenerating the 1.11 manifest leaves the shell on + whatever channel it was using. ### Fixed + - 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. + `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. - A remesh under `use_prior_polar` no longer resamples the refined sections' `SectionAero` surface tables down to whatever unrefined sections survive it. - `compute_refined_section_interpolation!` reblended contour, `cp` and `cf` from the - unrefined sections unconditionally while `aero_data` was preserved, so a wing rebuilt - onto fewer structural stations kept full-resolution polars but lost the surface tables - pressure integration reads. The reblend is now skipped when the polars are preserved - and the refined sections already carry tables. + `compute_refined_section_interpolation!` reblended contour, `cp` and `cf` from + the unrefined sections unconditionally while `aero_data` was preserved, so a + wing rebuilt onto fewer structural stations kept full-resolution polars but + lost the surface tables pressure integration reads. The reblend is now skipped + when the polars are preserved and the refined sections already carry tables. - The angle-of-attack correction (`correct_aoa=true`) no longer overwrites - `body_aero.AIC` with the aerodynamic-centre (LLT) matrix. It builds that matrix in - its own `AIC_aero_center` field, so `AIC` keeps holding the control-point matrix the - circulation was solved against, which is what anything reading it after a solve - expects. Costs a second `n_panels × n_panels × 3` buffer. -- `obj_to_yaml` no longer places sections on a wingtip that has closed to a point. - `station_indices` spreads its targets over the stations that still have a chord, - so the outermost section lands on the last sliceable one. A V3 mesh sliced with - the default `wingtip_distance` used to put a zero-chord section at each tip, - whose polar was `NaN` and took the whole solve with it; working around it meant - guessing a `wingtip_distance` large enough to skip past the tip. That workaround - is no longer the default: `wingtip_distance` is now `0.0`, an inset on top of the - trim for meshes whose slices just short of the tip are still too thin to analyse. -- Wing sections are normalized to `+y` to `-y` order on load (`normalize_span_order!`), - and by `refine!` for wings built through `add_section!`. Panel `y_airf` and `z_airf` - follow the order sections are stored in, so a geometry file written the other way - round inverted every panel normal, and a wing whose sections were replaced after its - panels were built (a structural remesh) inverted them mid-run, one panel at a time as - each crossed `spanwise_direction`. `obj_to_yaml` and `surfplan_to_aero_yaml` emit that + `body_aero.AIC` with the aerodynamic-centre (LLT) matrix. It builds that + matrix in its own `AIC_aero_center` field, so `AIC` keeps holding the + control-point matrix the circulation was solved against, which is what + anything reading it after a solve expects. Costs a second + `n_panels × n_panels × 3` buffer. +- `obj_to_yaml` no longer places sections on a wingtip that has closed to a + point. `station_indices` spreads its targets over the stations that still have + a chord, so the outermost section lands on the last sliceable one. A V3 mesh + sliced with the default `wingtip_distance` used to put a zero-chord section at + each tip, whose polar was `NaN` and took the whole solve with it; working + around it meant guessing a `wingtip_distance` large enough to skip past the + tip. That workaround is no longer the default: `wingtip_distance` is now + `0.0`, an inset on top of the trim for meshes whose slices just short of the + tip are still too thin to analyse. +- Wing sections are normalized to `+y` to `-y` order on load + (`normalize_span_order!`), and by `refine!` for wings built through + `add_section!`. Panel `y_airf` and `z_airf` follow the order sections are + stored in, so a geometry file written the other way round inverted every panel + normal, and a wing whose sections were replaced after its panels were built (a + structural remesh) inverted them mid-run, one panel at a time as each crossed + `spanwise_direction`. `obj_to_yaml` and `surfplan_to_aero_yaml` emit that order too; files of either order keep loading the same. -- Spanwise distribution plots put `+y` on the left, matching that order and the kite - seen from the front. -- The lofted airfoil skin in `plot_geometry` draws each section's contour on its own - panel edge (#256). The edge was picked from the wing's span order rather than from - the panel index, so a wing whose sections ran `-y` to `+y` got every rib drawn at its - neighbour's station, one tip bare and the other doubled, while the panels themselves - rendered correctly. +- Spanwise distribution plots put `+y` on the left, matching that order and the + kite seen from the front. +- The lofted airfoil skin in `plot_geometry` draws each section's contour on its + own panel edge (#256). The edge was picked from the wing's span order rather + than from the panel index, so a wing whose sections ran `-y` to `+y` got every + rib drawn at its neighbour's station, one tip bare and the other doubled, + while the panels themselves rendered correctly. ### 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. + +- `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. - `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. + `(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. - `is_show=true` draws into a window named after the plot title instead of into whichever window the backend last used, so a script showing several plots gets one window each and re-running it redraws them in place. `show_plot` takes the window `name` as a keyword. -- `examples/V3_kite.jl` builds its wing from `VSMSettings` and adds a second sweep on - polars generated from `V3_25.obj` with NeuralFoil (`NEURALFOIL`, on by default), so - the plots compare CAD-derived polars against the checked-in CFD tables. -- `bin/run_julia` starts Julia with `JULIA_NUM_THREADS=auto` unless the environment - already sets it. +- `examples/V3_kite.jl` builds its wing from `VSMSettings` and adds a second + sweep on polars generated from `V3_25.obj` with NeuralFoil (`NEURALFOIL`, on + by default), so the plots compare CAD-derived polars against the checked-in + CFD tables. +- `bin/run_julia` starts Julia with `JULIA_NUM_THREADS=auto` unless the + environment already sets it. ## VortexStepMethod v4.0.0 2026-08-03 ### Added -- `NeuralFoil`-based airfoil polar generation via the new `AirfoilAero` submodule - (`NeuralFoilSolver`, `XFoilSolver`, `analyze_section`, `analyze_sweep`, - `fit_kulfan_parameters`, `shrink_wrap`, `ShrinkWrap`) + +- `NeuralFoil`-based airfoil polar generation via the new `AirfoilAero` + submodule (`NeuralFoilSolver`, `XFoilSolver`, `analyze_section`, + `analyze_sweep`, `fit_kulfan_parameters`, `shrink_wrap`, `ShrinkWrap`) - `ObjAdapter` submodule: converts a 3D wing `.obj` mesh to the native YAML/CSV - geometry format (`obj_to_yaml`, `perpendicular_sections`, - `write_yaml`, `plot_slices_3d`, `plot_airfoils`) -- `SurfplanAdapter` submodule: `surfplan_to_aero_yaml` turns a SurfplanAdapter aero - export into the native pressure-ready YAML/CSV geometry (shared `generate_airfoils` - core with `ObjAdapter`) -- `ObjWing(obj_path[, dat_path]; Re, n_panels, aero_solver, remake, ...)` convenience - constructor restored for backward compatibility — internally calls + geometry format (`obj_to_yaml`, `perpendicular_sections`, `write_yaml`, + `plot_slices_3d`, `plot_airfoils`) +- `SurfplanAdapter` submodule: `surfplan_to_aero_yaml` turns a SurfplanAdapter + aero export into the native pressure-ready YAML/CSV geometry (shared + `generate_airfoils` core with `ObjAdapter`) +- `ObjWing(obj_path[, dat_path]; Re, n_panels, aero_solver, remake, ...)` + convenience constructor restored for backward compatibility — internally calls `ObjAdapter.obj_to_yaml` then `Wing`; `aero_solver` selects the polar backend (default `NeuralFoilSolver()`; pass `XFoilSolver()` for old behavior); - `remake=false` (default) reuses an existing `geometry.yaml` in `output_dir` to skip - expensive polar generation when only `n_panels` changes; set `remake=true` to force - regeneration -- Per-section surface aero tables: `SectionAero` (full contour + surface pressure - `cp` + skin friction `cf` per node over an `(α, δ)` grid), `section_surface`, - `read_section_aero` / `write_section_aero`; surface aero propagates through - `Section` and `Wing` and is spanwise-interpolated during `refine!` + `remake=false` (default) reuses an existing `geometry.yaml` in `output_dir` to + skip expensive polar generation when only `n_panels` changes; set + `remake=true` to force regeneration +- Per-section surface aero tables: `SectionAero` (full contour + surface + pressure `cp` + skin friction `cf` per node over an `(α, δ)` grid), + `section_surface`, `read_section_aero` / `write_section_aero`; surface aero + propagates through `Section` and `Wing` and is spanwise-interpolated during + `refine!` - `POLY` aero model for polynomial cl/cd/cm (exported) -- `lei_poly_coeffs(tube_diameter, camber)` (exported from `AirfoilAero`) returning - the Breukels α-polynomial cl/cd/cm coefficients +- `lei_poly_coeffs(tube_diameter, camber)` (exported from `AirfoilAero`) + returning the Breukels α-polynomial cl/cd/cm coefficients - Li/Gaunaa spanwise artificial viscosity (Li, Gaunaa, Pirrung & Lønbæk, TORQUE 2026) for the LOOP solver, stabilizing post-stall circulation distributions that otherwise develop non-physical sawtooth oscillations; opt-in via `is_with_artificial_viscosity` (default `false`) and `artificial_viscosity_factor` (default `0.035`) on the solver settings -- `crease_frac` wing setting (chordwise flap-hinge fraction, default `0.75`) - for drawing the δ-deflected plate/skin, readable from YAML +- `crease_frac` wing setting (chordwise flap-hinge fraction, default `0.75`) for + drawing the δ-deflected plate/skin, readable from YAML - `examples/V3_neuralfoil.jl` and `examples/obj_to_yaml_kite.jl` ### Changed + - BREAKING - `ObjWing` polar generation now uses NeuralFoil (via `ObjAdapter.obj_to_yaml`) instead of XFoil + a user-supplied `.dat` file. The constructor signature is backward-compatible (`dat_path` is accepted and @@ -227,26 +264,29 @@ `aero_solver=XFoilSolver()` to `ObjAdapter.obj_to_yaml` and call `Wing` directly. - OBJ-based wings are now built via `ObjAdapter.obj_to_yaml` + `Wing(yaml_path)` - instead of the old `ObjWing` pipeline (XFoil + single `.dat` file); - `ObjWing` is kept as a shim that accepts but ignores `dat_path` + instead of the old `ObjWing` pipeline (XFoil + single `.dat` file); `ObjWing` + is kept as a shim that accepts but ignores `dat_path` - BREAKING - `LEI_AIRFOIL_BREUKELS` is now a deprecated alias of `POLY`, and the built-in Breukels regression no longer runs at solve time. Sections must carry the α-polynomial coefficients directly instead of `(tube_diameter, camber)`; compute them up front with `lei_poly_coeffs(tube_diameter, camber)` - `Section` and `add_section!` accept an optional `section_aero` argument -- BREAKING - plotting is now Makie-only; the `VortexStepMethodMakieExt` extension - loads once a Makie backend and - [`MakieControlPlots`](https://github.com/OpenSourceAWE/MakieControlPlots.jl) are - available, so plotting code must now load `MakieControlPlots` (and drop - `set_plot_backend!`); `plot_section_polars` is rendered through `MakieControlPlots` +- BREAKING - plotting is now Makie-only; the `VortexStepMethodMakieExt` + extension loads once a Makie backend and + [`MakieControlPlots`](https://github.com/OpenSourceAWE/MakieControlPlots.jl) + are available, so plotting code must now load `MakieControlPlots` (and drop + `set_plot_backend!`); `plot_section_polars` is rendered through + `MakieControlPlots` ### Removed + - `ObjWing` as a standalone pipeline (replaced by `ObjAdapter`); the name is re-exported as a compatibility wrapper - `auto_rotation` helper (internal, removed from public API) -- `PanelGroupingMethod` enum (already removed in v3.0.0 — stale docs entry cleaned up) -- ControlPlots test run removed from CI (`plot-controlplots` arg) due to - a `libraqm`/HarfBuzz symbol conflict in the GitHub Actions environment +- `PanelGroupingMethod` enum (already removed in v3.0.0 — stale docs entry + cleaned up) +- ControlPlots test run removed from CI (`plot-controlplots` arg) due to a + `libraqm`/HarfBuzz symbol conflict in the GitHub Actions environment - BREAKING - the `ControlPlots` plotting extension, `examples_cp/`, and the `PythonCall`/Matplotlib setup (`bin/install_controlplots`, CondaPkg `LocalPreferences` defaults) @@ -258,6 +298,7 @@ ## VortexStepMethod v3.3.6 2026-06-13 ### Added + - `calc_forces!` and `solve_base!` (both exported): `solve!` is now `solve_base!` followed by `calc_forces!`, so a frozen circulation can be mapped to forces without re-running the nonlinear gamma solve (#245) @@ -266,10 +307,12 @@ kept as a thin wrapper (#246) ### Changed -- `calc_forces!` is now allocation-free in the per-step hot path - (preallocated `panel_area_dist` and `unrefined_count_dist` buffers) (#245) + +- `calc_forces!` is now allocation-free in the per-step hot path (preallocated + `panel_area_dist` and `unrefined_count_dist` buffers) (#245) ### Fixed + - 3D polar plotting (#245) - flaky Aqua `persistent_tasks` test now actually disabled via `persistent_tasks=false` (`()` did not disable it) (#246) @@ -277,19 +320,22 @@ ## VortexStepMethod v3.3.5 2026-06-05 ### Added + - `moment_coeff_unrefined_dist` field in `VSMSolution`: the summed `moment_frac`-referenced pitching-moment coefficient per unrefined section [-] ## VortexStepMethod v3.3.4 2026-05-31 ### Added -- `PlotBackend`, `MakieBackend`, `ControlPlotsBackend`, and - `set_plot_backend!` so applications can explicitly choose which plotting - extension the backend-agnostic plotting API should use -- `PythonCall` added as a weak dependency to support the `ControlPlots` - backend with PythonPlot + +- `PlotBackend`, `MakieBackend`, `ControlPlotsBackend`, and `set_plot_backend!` + so applications can explicitly choose which plotting extension the + backend-agnostic plotting API should use +- `PythonCall` added as a weak dependency to support the `ControlPlots` backend + with PythonPlot ### Changed + - backend-agnostic plotting wrappers now route through the active plotting backend, and each plotting extension initializes itself as the default only when no backend has been selected yet @@ -297,6 +343,7 @@ - improved `bin/install` and `bin/install_controlplots` scripts ### Fixed + - corrected projection onto core radius in `velocity_3D_bound_vortex!` and semi-infinite trailing vortex projection so that the radial direction is always measured from the filament axis, not from the origin (#241) @@ -307,6 +354,7 @@ ## VortexStepMethod v3.3.3 2026-05-21 ### Fixed + - `MakieExt` and `ControlPlotsExt` no longer both define `VortexStepMethod.plot_geometry` for the same type, resolving a method ambiguity when both extensions were loaded (#236) @@ -316,33 +364,39 @@ ## VortexStepMethod v3.3.2 2026-05-18 ### Changed + - use 2-arg version of atan to avoid possible NaN ## VortexStepMethod v3.3.1 2026-05-13 ### Changed + - `unrefined_deform!` linearly interpolates twist and TE deflection between unrefined sections and rotates each refined section about the average of its adjacent local airfoil normals (#234) ### Fixed -- `smooth_sqrt` in the solver hot loop keeps gradients defined at zero - velocity magnitude + +- `smooth_sqrt` in the solver hot loop keeps gradients defined at zero velocity + magnitude ## VortexStepMethod v3.3.0 2026-05-05 ### Added + - `ForwardDiff` compatibility, used by default in `linearize` (#232) - `backend` keyword argument for `linearize` - example `linearize_check.jl` comparing FiniteDiff and ForwardDiff tangents ### Changed + - core structs are parameterized on the scalar type `T` so dual numbers can propagate through them; public constructors are unchanged ## VortexStepMethod v3.2.0 2026-05-02 ### Added + - support for both CairoMakie and GLMakie in the example `menu()` via new `CairoMakie.activate()` / `GLMakie.activate()` entries (#222) - plots are saved as PDF when CairoMakie is active @@ -350,18 +404,21 @@ generated plots ### Changed + - examples write plots into the shared `output` folder instead of the working directory - example plot file names are sanitized: spaces replaced with `_` and `%` replaced with `pct` ### Fixed + - NONLIN solver no longer returns stale `gamma` on repeated `solve!` calls (#228) ## VortexStepMethod v3.1.3 2026-04-23 ### Fixed + - bug in `linearize` where `body_aero.va` was used instead of `body_aero._va`, causing incorrect initial velocity storage (#227) - `set_va!` no longer overwrites `_va` with a computed reference velocity when @@ -370,46 +427,63 @@ - linearize now correctly preserves and restores `omega` across perturbations ## VortexStepMethod v3.1.2 2025-04-20 + ### Added + - add back compat entry v2 for SciMLBase ## VortexStepMethod v3.1.1 2025-04-20 + ### Added + - add back compat entry v3 for RecursiveArrayTools ## VortexStepMethod v3.1.0 2025-04-19 ### Breaking + - `billowing_angle` replaced by `billowing_percentage` on `Wing` and `WingSettings` (percentage of arc length, not radians) - `billowing_angle_from_percentage()` removed - `BILLOWING` distribution now uses sinusoidal rotation instead of circular arc ### Added + - `billowing.jl` example comparing flat vs billowed V3 kite - Coarse V3 kite geometry, settings, and combined CFD polar data - `cl_over_cd` keyword for `plot_polars` and `plot_combined_analysis` -- the function `menu_cp()` can now be used to run the examples with the ControlPlots backend -- the script `bin/install` can and should be used to instantiate the project and all sub-projects after git checkout -- the files `Manifest-v1.11.toml.default` and `Manifest-v12.toml.default` for enhanced reproducibility -- the scripts `bin/jetls` and `bin/jetls_examples` for running a static check of the source code -- the script `test_bench.jl` for benchmarking the refinement method and measuring allocations +- the function `menu_cp()` can now be used to run the examples with the + ControlPlots backend +- the script `bin/install` can and should be used to instantiate the project and + all sub-projects after git checkout +- the files `Manifest-v1.11.toml.default` and `Manifest-v12.toml.default` for + enhanced reproducibility +- the scripts `bin/jetls` and `bin/jetls_examples` for running a static check of + the source code +- the script `test_bench.jl` for benchmarking the refinement method and + measuring allocations ### Changed + - Plot legends moved to shared horizontal legend at bottom of grid layouts - the script bin/run_julia now can be called with a script name as parameter -- fixed all JETLS warnings in the source code for improved performance and stability +- fixed all JETLS warnings in the source code for improved performance and + stability ### Fixed errors -- Domain error in elliptical gamma distribution when control points lie - outside the nominal span envelope + +- Domain error in elliptical gamma distribution when control points lie outside + the nominal span envelope ## VortexStepMethod v3.0.1 2025-04-04 ### Changed + - the file `CITATION.cff` - compat entry for RecursiveArrayTools + ### Added + - the file `.zenodo.json` ## VortexStepMethod v3.0.0 diff --git a/docs/src/functions.md b/docs/src/functions.md index 94b45af4..dc7fdcf2 100644 --- a/docs/src/functions.md +++ b/docs/src/functions.md @@ -47,7 +47,13 @@ chord_line control_point_deflection panel_kulfan_parameters refresh_live_polars! +deform_live_shapes! +apply_live_shapes! polar_drift +refresh_live_pressure! +live_surface_friction! +contour_shape_matrix +live_shape_offset! generate_aero_matrices generate_polar_from_coordinates generate_polar_from_dat diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 0691765a..e657184c 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -98,8 +98,10 @@ update_non_deformed_sections! ### Aerodynamic data and Cp ```@docs calculate_new_aero_data -set_taylor_polar! -taylor_value +set_sampled_polar! +decode_surface_pressure +contour_pressure +sampled_value assemble_polar_matrix load_matrix_polar_data read_aero_matrix diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 57b0ebb9..62b7a859 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -36,7 +36,7 @@ export calculate_projected_area, calculate_span export MVec3 export LLT, Model, VSM -export AeroModel, INVISCID, POLY, LEI_AIRFOIL_BREUKELS, POLAR_MATRICES, POLAR_VECTORS, TAYLOR +export AeroModel, INVISCID, POLY, LEI_AIRFOIL_BREUKELS, POLAR_MATRICES, POLAR_VECTORS, SAMPLED export KulfanParameters export BILLOWING, COSINE, LINEAR, PanelDistribution, SPLIT_PROVIDED, UNCHANGED export ELLIPTIC, InitialGammaDistribution, ZEROS @@ -238,7 +238,7 @@ Enumeration of the implemented wing types. @enum WingType RECTANGULAR CURVED ELLIPTICAL """ - AeroModel `POLY` `POLAR_VECTORS` `POLAR_MATRICES` `INVISCID` `TAYLOR` + AeroModel `POLY` `POLAR_VECTORS` `POLAR_MATRICES` `INVISCID` `SAMPLED` Enumeration of the implemented aerodynamic models. See also: [AeroData](@ref) @@ -248,10 +248,11 @@ Enumeration of the implemented aerodynamic models. See also: [AeroData](@ref) - `POLAR_VECTORS`: Polar vectors as function of alpha (lookup tables with interpolation) - `POLAR_MATRICES`: Polar matrices as function of alpha and delta (lookup tables with interpolation) - INVISCID -- `TAYLOR`: order-N polynomial in `α - α_ref` [rad] per panel, the local expansion a - live polar source refits each solve, see [`refresh_live_polars!`](@ref - VortexStepMethod.AirfoilAero.refresh_live_polars!). Valid only inside its fit - window, so it carries no stall model and ignores `delta`. +- `SAMPLED`: cl/cd/cm sampled at ascending angles of attack per panel, interpolated + between them and held flat past either end. This is what a live polar source writes + each solve, see [`refresh_live_polars!`](@ref + VortexStepMethod.AirfoilAero.refresh_live_polars!). Samples are the polar, so a stall + knee inside the sampled range is represented rather than smoothed. Ignores `delta`. `LEI_AIRFOIL_BREUKELS` is a deprecated alias of `POLY`. @@ -262,7 +263,7 @@ where `alpha` is the angle of attack, `delta` is trailing edge angle. POLAR_VECTORS POLAR_MATRICES INVISCID - TAYLOR + SAMPLED end """ @@ -336,8 +337,7 @@ abstract type AbstractWing{T} end Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}}, - Tuple{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, - Tuple{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}} + Tuple{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}} } Union of different definitions of the aerodynamic properties of a wing section. See also: [AeroModel](@ref) @@ -345,8 +345,9 @@ Union of different definitions of the aerodynamic properties of a wing section. - (`cl_coeffs`, `cd_coeffs`, `cm_coeffs`) α-polynomial coefficients for `POLY` - (`alpha_range`, `cl_vector`, `cd_vector`, `cm_vector`) for `POLAR_VECTORS` - (`alpha_range`, `delta_range`, `cl_matrix`, `cd_matrix`, `cm_matrix`) for `POLAR_MATRICES` - - (`alpha_ref`, `cl_coeffs`, `cd_coeffs`, `cm_coeffs`) for `TAYLOR`, the expansion - point [rad] and the ascending coefficients of the polynomial in `α - alpha_ref` + +`SAMPLED` carries no section-level data: it is written onto a panel at run time by a +live polar source, never read from a section. where `alpha` is the angle of attack [rad], `delta` is trailing edge angle [rad], `cl` the lift coefficient, `cd` the drag coefficient and `cm` the pitching moment coefficient. The camber of a kite refers to @@ -358,8 +359,7 @@ const AeroData = Union{ Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}}, - Tuple{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, - Tuple{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}} + Tuple{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}} } const PACKAGE_ROOT = normpath(joinpath(@__DIR__, "..")) diff --git a/src/airfoil_aero/AirfoilAero.jl b/src/airfoil_aero/AirfoilAero.jl index abf345dc..d02b1acd 100644 --- a/src/airfoil_aero/AirfoilAero.jl +++ b/src/airfoil_aero/AirfoilAero.jl @@ -8,7 +8,7 @@ using NPZ using Xfoil using Printf: @sprintf using ..VortexStepMethod: SectionAero, interpolate_matrix_nans!, delta_suffix, - write_node_rows, section_surface, set_taylor_polar!, + write_node_rows, section_surface, set_sampled_polar!, KulfanParameters include("kulfan.jl") @@ -35,7 +35,9 @@ export neuralfoil_aero, neuralfoil_section export KulfanBasis, deform_kulfan, control_point_deflection export chord_residual, chord_line export LivePolarSettings, LivePolars, panel_kulfan_parameters -export refresh_live_polars!, polar_drift +export refresh_live_polars!, refresh_live_pressure!, live_surface_friction! +export contour_shape_matrix, live_shape_offset! +export polar_drift export AbstractAirfoilSolver, XFoilSolver, NeuralFoilSolver export SectionSolution, DeformedSection, deform_section, analyze_section, analyze_sweep export create_2d_polars, generate_aero_matrices, generate_section_aero, lei_poly_coeffs diff --git a/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl b/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl index 9a625600..6a807d39 100644 --- a/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl +++ b/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl @@ -47,9 +47,8 @@ Assemble a full-contour [`SectionSolution`](@ref) for case `i` of a NeuralFoil s `(x, y)` split at the leading edge `le`, carrying the precomputed `cf`. """ function neuralfoil_contour_solution(alpha, res, i, x, y, le, cf) - up = linear_interpolation(res.x, res.cp_upper[:, i]; extrapolation_bc=Line()) - lo = linear_interpolation(res.x, res.cp_lower[:, i]; extrapolation_bc=Line()) - cp = [(k <= le ? up : lo)(clamp(x[k], 0.0, 1.0)) for k in eachindex(x)] + cp = contour_pressure(res.x, view(res.cp_upper, :, i), + view(res.cp_lower, :, i), x, le) return SectionSolution(alpha, res.cl[i], res.cd[i], res.cm[i], res.confidence[i], x, y, cp, cf) end diff --git a/src/airfoil_aero/deform.jl b/src/airfoil_aero/deform.jl index 08b61779..c0f551fa 100644 --- a/src/airfoil_aero/deform.jl +++ b/src/airfoil_aero/deform.jl @@ -1,9 +1,17 @@ """ - KulfanBasis(; n_stations=60, n_weights=8) + KulfanBasis(; n_stations=60, n_weights=8, ridge=1e-3) The fixed CST basis a shape deformation is projected onto: chord stations `x` in -`[0, 1]` and the pseudoinverse `projection` of the class-times-Bernstein matrix -`C(x)·B(x)` those stations span. +`[0, 1]` and the ridge-regularised inverse `projection` of the class-times-Bernstein +matrix `C(x)·B(x)` those stations span. + +`ridge` is the Tikhonov weight, relative to the basis' own largest singular value, that +keeps the projection from answering a deflection it cannot represent with weights far +larger than the airfoil they correct. A plain pseudoinverse has no such bound: a kinked +deflection — a strut buckling is one — lands on the basis' weakest directions and comes +back amplified a hundredfold and alternating in sign, which is not an airfoil. The ridge +trades a small, measured under-response on deflections the basis *can* hold for a bounded +answer on the ones it cannot. CST is linear in its weights, so a surface displacement is a matvec against this constant matrix — never a refit. Refitting inside a loop is not an option: the Kulfan @@ -18,19 +26,25 @@ the representable residual. struct KulfanBasis "Chord stations the deflection is sampled on, ascending in `[0, 1]`." x::Vector{Float64} - "`n_weights × length(x)` pseudoinverse of `C(x)·B(x)`, mapping deflection to weights." + "`n_weights × length(x)` regularised inverse of `C(x)·B(x)`, mapping deflection to weights." projection::Matrix{Float64} "Number of CST weights per surface, matching the airfoil being deformed." n_weights::Int + "Tikhonov weight the projection was built with, relative to the basis' largest singular value." + ridge::Float64 end -function KulfanBasis(; n_stations::Int=60, n_weights::Int=8) +function KulfanBasis(; n_stations::Int=60, n_weights::Int=8, ridge::Real=1e-3) n_stations > n_weights || throw(ArgumentError( "KulfanBasis needs more stations than weights; got $n_stations and $n_weights.")) + ridge >= 0 || throw(ArgumentError("KulfanBasis ridge must not be negative.")) theta = range(0, pi, n_stations) x = @. (1 - cos(theta)) / 2 shape = class_function(x) .* bernstein_basis(x, n_weights - 1) - return KulfanBasis(collect(x), pinv(shape), n_weights) + scale = maximum(svdvals(shape))^2 + projection = iszero(ridge) ? pinv(shape) : + (shape' * shape + ridge * scale * I) \ shape' + return KulfanBasis(collect(x), projection, n_weights, Float64(ridge)) end """ diff --git a/src/airfoil_aero/live_polar.jl b/src/airfoil_aero/live_polar.jl index 7edd7bae..65e37515 100644 --- a/src/airfoil_aero/live_polar.jl +++ b/src/airfoil_aero/live_polar.jl @@ -1,30 +1,25 @@ """ - LivePolarSettings(; order=2, half_window=deg2rad(4), n_samples=5, - model_size="xlarge", weights_dir=nothing, n_crit=9.0) + LivePolarSettings(; offsets=deg2rad.(-12.0:3.0:12.0), model_size="xlarge", + weights_dir=nothing, n_crit=9.0) -How a live polar is sampled and fitted. Every solve, each panel's deformed shape is -evaluated at `n_samples` angles of attack spanning `α_ref ± half_window` and a -polynomial of `order` is least-squares fitted through them, becoming the panel's -`TAYLOR` polar. +How a live polar is sampled. Every solve, each panel's deformed shape is evaluated at +its reference angle of attack plus each of `offsets`, and those values become the +panel's `SAMPLED` polar directly — there is no fit in between. -Order and window belong together. Order 2 over `±4°` tracks a full-polar solve to -within a few tenths of a percent and bends with the stall knee where a straight line -cuts across it; higher orders carry large coefficients that diverge as soon as the -solve steps outside the window, which is why the window is a correctness limit and not -an accuracy knob — see [`polar_drift`](@ref). +Sampling rather than fitting is what lets a panel hold a stall. A polynomial over the +same window averages the knee into a slope and past the peak returns a lift slope of +the wrong sign, which is a divergence the solve cannot recover from; sampled values +are the polar at their own angles and reproduce whatever shape lies between them. -Fitting is by least squares over the window, never by finite differences at a tabulated -source's own knot spacing: a piecewise-linear polar has zero curvature inside a segment -and a spike at every knot, so a difference at the knot spacing returns pure -discretization artefact. +The offsets set both what the polar resolves and how far the solve may move before it +reads a flat end (see [`polar_drift`](@ref)): their spacing is the resolution a stall +knee is caught at, their reach is the room the solve has. The default spacing matches +the 3° grid the offline polar tables are generated on, over a range wide enough that a +tip crossing into stall stays inside it. """ @with_kw struct LivePolarSettings - "Polynomial order in `α - α_ref`; 2 unless the window is narrowed to match." - order::Int = 2 - "Half width [rad] of the fit window, also the drift the solve may not leave." - half_window::Float64 = deg2rad(4.0) - "Angles of attack sampled per panel per refresh." - n_samples::Int = 5 + "Angles of attack [rad] off the reference angle, ascending and straddling zero." + offsets::Vector{Float64} = deg2rad.(collect(-12.0:3.0:12.0)) "NeuralFoil network size." model_size::String = "xlarge" "Directory holding the network weights; `nothing` takes the packaged ones." @@ -37,13 +32,13 @@ end LivePolars(base; settings=LivePolarSettings(), n_stations=60) Live polar source for a wing whose panels each carry one undeformed airfoil in `base`. -Holds the fixed CST basis, the per-panel expansion points and the network input scratch. +Holds the fixed CST basis, the per-panel reference angles and the network input scratch. A refresh is dominated by the forward pass itself; the deformation is a matvec against the constant basis, about half a microsecond a panel. Drive it with [`refresh_live_polars!`](@ref). """ mutable struct LivePolars - "Sampling and fitting configuration." + "Sampling configuration." settings::LivePolarSettings "The fixed CST basis every panel's deflection is projected onto." basis::KulfanBasis @@ -51,31 +46,34 @@ mutable struct LivePolars base::Vector{KulfanParameters} "Deformed Kulfan parameters per panel, rewritten every refresh." deformed::Vector{KulfanParameters} - "Expansion point [rad] the last fit was built about, per panel." + "Reference angle [rad] the last refresh sampled about, per panel." alpha_ref::Vector{Float64} - "Sample offsets [rad] off the expansion point, shared by every panel." - offsets::Vector{Float64} - "`(order + 1) × n_samples` least-squares solver for the shared offset grid." - fit::Matrix{Float64} + "One panel's sample angles [rad], rebuilt per panel inside a refresh." + knots::Vector{Float64} "`25 × (n_panels · n_samples)` network input scratch." inputs::Matrix{Float32} + "`25 × n_panels` network input scratch for the surface-pressure pass." + pressure_inputs::Matrix{Float32} end function LivePolars(base::AbstractVector{KulfanParameters}; settings::LivePolarSettings=LivePolarSettings(), n_stations::Int=60) - settings.n_samples > settings.order || throw(ArgumentError( - "LivePolars needs more samples than the fit order; got " * - "$(settings.n_samples) and $(settings.order).")) - offsets = collect(range(-settings.half_window, settings.half_window, - settings.n_samples)) - vandermonde = [d^(k - 1) for d in offsets, k in 1:(settings.order + 1)] + offsets = settings.offsets + length(offsets) >= 2 || throw(ArgumentError( + "LivePolars needs at least two sample offsets; got $(length(offsets)).")) + issorted(offsets) || throw(ArgumentError( + "LivePolars needs ascending sample offsets.")) + first(offsets) <= 0 <= last(offsets) || throw(ArgumentError( + "LivePolars sample offsets must straddle zero, so the reference angle lies " * + "inside the sampled range; got $(rad2deg.(extrema(offsets))) deg.")) n_panels = length(base) return LivePolars(settings, KulfanBasis(; n_stations, n_weights=length(base[1].upper_weights)), - collect(base), collect(base), zeros(n_panels), offsets, - pinv(vandermonde), - zeros(Float32, 25, n_panels * settings.n_samples)) + collect(base), collect(base), zeros(n_panels), + zeros(length(offsets)), + zeros(Float32, 25, n_panels * length(offsets)), + zeros(Float32, 25, n_panels)) end """ @@ -103,26 +101,67 @@ end """ polar_drift(live, alpha) -> Float64 -How far the largest panel angle of attack has drifted off the expansion point its -polar was fitted about, as a fraction of the fit window. Above `1` a panel is being -evaluated outside the window and the fit must be rebuilt before its answer is used, so -this is the guard [`refresh_live_polars!`](@ref) leaves to the caller's solve loop. +How far the largest panel angle of attack has drifted off the reference angle its polar +was sampled about, as a fraction of the reach the samples give it. Above `1` a panel is +being evaluated past the last sample, where the polar is held flat and says nothing +about how the panel is really behaving, so this is the diagnostic a caller's solve loop +reports or refreshes on. Asymmetric offsets are measured by their shorter side. """ function polar_drift(live::LivePolars, alpha::AbstractVector) length(alpha) == length(live.alpha_ref) || throw(ArgumentError( "polar_drift: $(length(alpha)) angles for $(length(live.alpha_ref)) panels.")) - return maximum(abs.(alpha .- live.alpha_ref)) / live.settings.half_window + offsets = live.settings.offsets + reach = min(-first(offsets), last(offsets)) + return maximum(abs.(alpha .- live.alpha_ref)) / reach +end + +""" + deform_live_shapes!(live::LivePolars, deflection) -> Vector{KulfanParameters} + +Deform every base airfoil by its camber increment and store the result in +`live.deformed`, which is returned. `nothing` leaves the base shapes in place. +The deformation is an analytic perturbation of one fixed weight vector, so it +never refits and never inherits the non-uniqueness of a fit. +""" +function deform_live_shapes!(live::LivePolars, deflection) + for i in eachindex(live.base) + live.deformed[i] = isnothing(deflection) ? live.base[i] : + deform_kulfan(live.basis, live.base[i], deflection[i]) + end + return live.deformed +end + +""" + apply_live_shapes!(live::LivePolars, panels; deflection=nothing) + +Write each panel's deformed airfoil into `live_shape` without touching its polar, +so the shape a panel reports is the one its current deformation gives. Replaying a +log re-derives the drawn airfoil this way: the frame is never solved, so its polars +would say nothing, and the network pass they cost is the whole price of a refresh. +""" +function apply_live_shapes!(live::LivePolars, panels; deflection=nothing) + length(panels) == length(live.base) || throw(ArgumentError( + "apply_live_shapes!: $(length(panels)) panels for $(length(live.base)) " * + "base airfoils.")) + deform_live_shapes!(live, deflection) + for (panel, shape) in zip(panels, live.deformed) + panel.live_shape = shape + end + return nothing end """ refresh_live_polars!(live, panels, alpha_ref, reynolds; deflection=nothing) -> Float64 -Refit every panel's polar from its current shape and write it in as a `TAYLOR` polar -(see [`set_taylor_polar!`](@ref VortexStepMethod.set_taylor_polar!)). Per panel: -deform the base airfoil by `deflection` (a chord-normalized deflection on -`live.basis.x`, or `nothing` to keep the base shape), evaluate NeuralFoil at -`alpha_ref ± half_window`, and least-squares fit the polynomial. +Regenerate every panel's polar from its current shape and write it in as a `SAMPLED` +polar (see [`set_sampled_polar!`](@ref VortexStepMethod.set_sampled_polar!)). Per panel: +deform the base airfoil by `deflection` (a chord-normalized deflection on `live.basis.x`, +or `nothing` to keep the base shape), evaluate NeuralFoil at `alpha_ref .+ offsets`, and +hand those values straight to the panel. + +The write is in place — same knot count, same vectors — so a refresh every solve costs +the forward pass and nothing else. Each panel keeps the deformed shape it was evaluated at as its `live_shape`, so a plot draws the airfoil the network actually saw. @@ -133,7 +172,7 @@ the batch, a value near zero meaning the deformed shape has left the region the was trained on. Call it after the mesh has been rebuilt from the structure and before the solve: a mesh -rebuild re-seeds each panel's aero from its section, which would drop the fit. +rebuild re-seeds each panel's aero from its section, which would drop the polar. """ function refresh_live_polars!(live::LivePolars, panels, alpha_ref, reynolds; deflection=nothing) @@ -144,13 +183,13 @@ function refresh_live_polars!(live::LivePolars, panels, alpha_ref, reynolds; alpha_vec, re_vec = per_panel(alpha_ref), per_panel(reynolds) live.alpha_ref .= alpha_vec - n_samples = live.settings.n_samples + offsets = live.settings.offsets + n_samples = length(offsets) + deform_live_shapes!(live, deflection) for i in 1:n_panels - live.deformed[i] = isnothing(deflection) ? live.base[i] : - deform_kulfan(live.basis, live.base[i], deflection[i]) for k in 1:n_samples fill_case_input!(live.inputs, (i - 1) * n_samples + k, live.deformed[i], - rad2deg(alpha_vec[i] + live.offsets[k]), re_vec[i], + rad2deg(alpha_vec[i] + offsets[k]), re_vec[i], live.settings.n_crit, 1.0, 1.0) end end @@ -161,9 +200,114 @@ function refresh_live_polars!(live::LivePolars, panels, alpha_ref, reynolds; for i in 1:n_panels samples = ((i - 1) * n_samples + 1):(i * n_samples) - set_taylor_polar!(panels[i], alpha_vec[i], live.fit * cl[samples], - live.fit * cd[samples], live.fit * cm[samples]; - window=live.settings.half_window, shape=live.deformed[i]) + live.knots .= alpha_vec[i] .+ offsets + set_sampled_polar!(panels[i], live.knots, view(cl, samples), + view(cd, samples), view(cm, samples); + shape=live.deformed[i]) end return minimum(confidence) end + +""" + refresh_live_pressure!(cp, live, contour_x, leading_edge, alpha, reynolds) + -> Vector{Vector{Float64}} + +Fill `cp[i]` with the surface pressure of panel `i`'s current deformed shape at +`alpha[i]`, resampled onto the contour nodes `contour_x[i]` (see +[`contour_pressure`](@ref)). `cp` is written in place and returned. + +This is the pressure half of a live polar. [`refresh_live_polars!`](@ref) makes the +panel forces follow the deformed shape; without this the pattern that spreads those +forces over the structure would still come from the undeformed section, so a +deformation would change how hard a panel pulls but not where it pulls. + +One batched forward pass over all panels, at the converged angle of attack rather +than at the sampled ones — a panel's `Cp` is wanted at exactly one angle, and +evaluating there is both cheaper than storing the samples and exact. Call it after +the solve has converged, with `contour_x` and `leading_edge` from the panel contours +the traction pattern is indexed on. Deform the shapes first, which +[`refresh_live_polars!`](@ref) already did for this solve. +""" +function refresh_live_pressure!(cp, live::LivePolars, contour_x, leading_edge, + alpha, reynolds) + n_panels = length(live.base) + length(cp) == length(contour_x) == length(leading_edge) == n_panels || + throw(ArgumentError("refresh_live_pressure!: $(length(cp)) pressure and " * + "$(length(contour_x)) contour entries for $n_panels panels.")) + per_panel(v) = v isa Number ? fill(float(v), n_panels) : collect(float.(v)) + alpha_vec, re_vec = per_panel(alpha), per_panel(reynolds) + for i in 1:n_panels + fill_case_input!(live.pressure_inputs, i, live.deformed[i], + rad2deg(alpha_vec[i]), re_vec[i], live.settings.n_crit, + 1.0, 1.0) + end + model = load_neuralfoil_model(live.settings.model_size; + weights_dir=live.settings.weights_dir) + station_x, cp_upper, cp_lower = decode_surface_pressure( + fused_output(live.pressure_inputs, model)) + for i in 1:n_panels + cp[i] .= contour_pressure(station_x, view(cp_upper, :, i), + view(cp_lower, :, i), contour_x[i], + leading_edge[i]) + end + return cp +end + +""" + live_surface_friction!(cf, contour_x, reynolds) -> Vector{Vector{Float64}} + +Fill `cf[i]` with the skin friction of panel `i`'s contour nodes at `reynolds[i]`, +by the same flat-plate closure ([`flat_plate_cf`](@ref)) the offline tables carry — +NeuralFoil does not predict skin friction. It depends on chord fraction and Reynolds +only, not on the shape or the angle of attack, so the live value differs from the +tabulated one purely by being at the panel's own flight Reynolds instead of the one +the tables were generated at. +""" +function live_surface_friction!(cf, contour_x, reynolds) + length(cf) == length(contour_x) || throw(ArgumentError( + "live_surface_friction!: $(length(cf)) friction and " * + "$(length(contour_x)) contour entries.")) + re_vec = reynolds isa Number ? fill(float(reynolds), length(cf)) : + collect(float.(reynolds)) + for i in eachindex(cf) + cf[i] .= (flat_plate_cf(clamp(x, 0.0, 1.0), re_vec[i]) + for x in contour_x[i]) + end + return cf +end + +""" + contour_shape_matrix(contour_x, n_weights) -> Matrix{Float64} + +The CST shape matrix `C(x)·B(x)` at chord fractions `contour_x`, mapping a change in +Kulfan weights to the normal offset it produces there. Build it once per contour: the +chord fractions of a panel's contour nodes never move, only the weights do. +""" +function contour_shape_matrix(contour_x, n_weights) + x = clamp.(collect(float.(contour_x)), 0.0, 1.0) + return class_function(x) .* bernstein_basis(x, n_weights - 1) +end + +""" + live_shape_offset!(offset, live::LivePolars, shape) -> Vector{Vector{Float64}} + +Fill `offset[i]` with the normal offset, over chord, that panel `i`'s current +deformation adds to its contour, given the panel's `shape` matrix from +[`contour_shape_matrix`](@ref). Written in place and returned. + +A deflection deforms both surfaces by the same camber increment, so one offset serves +the upper and lower halves of a contour alike and the surface split +[`contour_pressure`](@ref) needs does not arise here. Added to a reference contour it +gives the surface the network was actually evaluated on, which is what a traction +pattern has to be draped over for its normals and segment areas to mean anything. +""" +function live_shape_offset!(offset, live::LivePolars, shape) + length(offset) == length(shape) == length(live.base) || throw(ArgumentError( + "live_shape_offset!: $(length(offset)) offset and $(length(shape)) shape " * + "entries for $(length(live.base)) panels.")) + for i in eachindex(live.base) + mul!(offset[i], shape[i], + live.deformed[i].upper_weights .- live.base[i].upper_weights) + end + return offset +end diff --git a/src/airfoil_aero/neuralfoil.jl b/src/airfoil_aero/neuralfoil.jl index ea028f57..afd9928e 100644 --- a/src/airfoil_aero/neuralfoil.jl +++ b/src/airfoil_aero/neuralfoil.jl @@ -407,18 +407,46 @@ function neuralfoil_section(params::KulfanParameters, alpha, Re; n_crit=9.0, xtr_upper=1.0, xtr_lower=1.0) y = neuralfoil_fused_output(params, alpha, Re; model_size, weights_dir, n_crit, xtr_upper, xtr_lower) - N = (size(y, 1) - 6) ÷ 6 - upper_ue = y[(7 + 2N):(6 + 3N), :] - lower_ue = y[(7 + 5N):(6 + 6N), :] alpha_vec = alpha isa Number ? [Float64(alpha)] : Float64.(collect(alpha)) + station_x, cp_upper, cp_lower = decode_surface_pressure(y) return (; alpha = alpha_vec, cl = Vector{Float64}(y[2, :] ./ 2), cd = Vector{Float64}(clamp.(exp.((y[3, :] .- 2) .* 2), 0.0, 1.0)), cm = Vector{Float64}(y[4, :] ./ 20), confidence = Vector{Float64}(sigmoid.(y[1, :])), - x = compute_optimal_x_points(N), - cp_upper = Matrix{Float64}(1 .- upper_ue .^ 2), - cp_lower = Matrix{Float64}(1 .- lower_ue .^ 2)) + x = station_x, cp_upper, cp_lower) +end + +""" + decode_surface_pressure(y) -> (station_x, cp_upper, cp_lower) + +Surface pressure of a fused output matrix, reconstructed from the predicted +edge-velocity ratios as `Cp = 1 - (ue/vinf)^2`. `station_x` is NeuralFoil's own +`N` fixed stations, and each `Cp` matrix is `N × n_cases`. +""" +function decode_surface_pressure(y::AbstractMatrix) + N = (size(y, 1) - 6) ÷ 6 + upper_ue = y[(7 + 2N):(6 + 3N), :] + lower_ue = y[(7 + 5N):(6 + 6N), :] + return (compute_optimal_x_points(N), + Matrix{Float64}(1 .- upper_ue .^ 2), + Matrix{Float64}(1 .- lower_ue .^ 2)) +end + +""" + contour_pressure(station_x, cp_upper, cp_lower, contour_x, leading_edge) + -> Vector{Float64} + +`Cp` at every node of a closed contour, taken from the two surface distributions +NeuralFoil predicts at `station_x`. Nodes up to `leading_edge` read the upper +surface and the rest the lower, which is the Selig ordering the contour is stored +in — trailing edge, over the top to the nose, back along the bottom. +""" +function contour_pressure(station_x, cp_upper, cp_lower, contour_x, leading_edge) + upper = linear_interpolation(station_x, cp_upper; extrapolation_bc=Line()) + lower = linear_interpolation(station_x, cp_lower; extrapolation_bc=Line()) + return [(k <= leading_edge ? upper : lower)(clamp(contour_x[k], 0.0, 1.0)) + for k in eachindex(contour_x)] end """ diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index f8ee8f7a..a97b007e 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -224,9 +224,9 @@ function calculate_stall_angle_list!(stall_angles::AbstractVector, # Default stall angle if none found panel_stall = stall_angle_if_none_detected - # A local expansion says nothing outside its fit window, so its curvature must - # not be read as a stall peak far from the operating point. - if panel.aero_model == TAYLOR + # A live polar only spans the window it was sampled over, so neither its + # curvature nor its flat ends may be read as a stall peak outside that. + if panel.aero_model == SAMPLED stall_angles[idx] = panel_stall continue end diff --git a/src/panel.jl b/src/panel.jl index e29d1f73..0266eab0 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -36,8 +36,9 @@ Represents a panel in a vortex step method simulation. All points and vectors ar ): Panel filaments, see: [BoundFilament](@ref) - `delta`::T=0: flap trailing-edge deflection [rad] - `crease_frac`::T=0: chordwise flap-hinge fraction (0–1); 0 disables the plate kink -- `alpha_ref`::Float64=0: expansion point [rad] of the `TAYLOR` coefficients -- `alpha_window`::Float64=0: half width [rad] the `TAYLOR` fit is valid over; 0 = unbounded +- `alpha_ref`::Float64=0: reference angle [rad] of the `SAMPLED` polar +- `alpha_window`::Float64=0: half width [rad] that polar reaches; 0 = unbounded +- `alpha_knots`::Vector{Float64}=Float64[]: ascending angles [rad] a `SAMPLED` polar holds values at - `live_shape`::Union{Nothing, KulfanParameters}=nothing: the deformed airfoil the polar was generated from """ @with_kw mutable struct Panel{T, CL, CD, CM, SA} @@ -75,6 +76,7 @@ Represents a panel in a vortex step method simulation. All points and vectors ar crease_frac::T = zero(T) alpha_ref::Float64 = 0.0 alpha_window::Float64 = 0.0 + alpha_knots::Vector{Float64} = Float64[] live_shape::Union{Nothing, KulfanParameters} = nothing end @@ -222,16 +224,6 @@ function init_aero!(panel::Panel, section_1::Section, section_2::Section; panel.cl_coeffs = (c1[1] .+ c2[1]) ./ 2 panel.cd_coeffs = (c1[2] .+ c2[2]) ./ 2 panel.cm_coeffs = (c1[3] .+ c2[3]) ./ 2 - elseif panel.aero_model == TAYLOR - c1, c2 = section_1.aero_data, section_2.aero_data - (c1 isa Tuple{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}} && - c2 isa Tuple{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}}) || - throw(ArgumentError("TAYLOR requires aero_data = " * - "(alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs).")) - all(length.(c1[2:4]) .== length.(c2[2:4])) || - throw(ArgumentError("TAYLOR coefficient vectors must have equal length.")) - set_taylor_polar!(panel, (c1[1] + c2[1]) / 2, (c1[2] .+ c2[2]) ./ 2, - (c1[3] .+ c2[3]) ./ 2, (c1[4] .+ c2[4]) ./ 2) elseif !(panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES, INVISCID)) throw(ArgumentError("Unsupported aero model: $(panel.aero_model)")) end @@ -241,61 +233,69 @@ function init_aero!(panel::Panel, section_1::Section, section_2::Section; end """ - taylor_value(coeffs, delta_alpha, window) + sampled_value(knots, values, alpha) -A `TAYLOR` polar's coefficient at `delta_alpha = α - α_ref` [rad], continued linearly -beyond `±window`. A local expansion's arms diverge fast outside the range it was fitted -over — an order-2 fit is worse than useless a few degrees out — so past the edge the -polynomial's own value and slope there carry it on instead. That keeps a solve stepping -outside the window convergent and honest about being extrapolated, rather than chasing a -parabola to infinity; the caller's drift guard is what puts the fit back where the solve -went. `window = 0` leaves the polynomial unbounded. +A `SAMPLED` polar's coefficient at `alpha` [rad]: linear between the knots it was +sampled on, and the end value beyond either end. Flat ends are the point — the +polar of a stalled section has a knee that no polynomial holds, and a coefficient +carried on past the samples must stay bounded rather than follow a slope out. -The continuation is a straight line, so it says nothing about stall and can carry a -coefficient somewhere physically impossible. [`calculate_cd`](@ref) therefore floors -its result at zero: every other coefficient may be extrapolated, but a negative drag -would feed energy into whatever reads it. +`knots` is ascending; a single knot is a constant polar. """ -@inline function taylor_value(coeffs, delta_alpha, window) - if window > 0 && abs(delta_alpha) > window - edge = delta_alpha >= 0 ? window : -window - slope = sum((k - 1) * coeffs[k] * edge^(k - 2) for k in 2:length(coeffs); - init = 0.0) - return evalpoly(edge, coeffs) + slope * (delta_alpha - edge) +@inline function sampled_value(knots, values, alpha) + n = length(knots) + n == 0 && throw(ArgumentError("A SAMPLED polar has no knots.")) + alpha <= knots[1] && return values[1] + alpha >= knots[n] && return values[n] + hi = 2 + @inbounds while alpha > knots[hi] + hi += 1 + end + @inbounds begin + span = knots[hi] - knots[hi - 1] + blend = span > 0 ? (alpha - knots[hi - 1]) / span : zero(alpha) + return values[hi - 1] + blend * (values[hi] - values[hi - 1]) end - return evalpoly(delta_alpha, coeffs) end """ - set_taylor_polar!(panel, alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs; window=0.0, - shape=nothing) + set_sampled_polar!(panel, alphas, cl, cd, cm; shape=nothing) + +Overwrite a panel's local polar with values sampled at `alphas` [rad], ascending. +The panel's aero model is set to `SAMPLED`. + +Written in place: the knots and the three value vectors reuse the panel's own +`alpha_knots` and `cl_coeffs`/`cd_coeffs`/`cm_coeffs` storage whenever the sample +count is unchanged, which is what lets a live polar source refresh every solve +without allocating. The vectors hold sampled values here rather than polynomial +coefficients — the aero model is what says which. -Overwrite a `TAYLOR` panel's local polar: the expansion point `alpha_ref` [rad], the -ascending coefficients of the polynomials in `α - alpha_ref`, and the half width -`window` [rad] the fit is valid over, past which it is continued linearly (see -[`taylor_value`](@ref)). Copies into the panel's existing coefficient vectors when the -order is unchanged, so a live polar source can refit every solve without allocating. -The panel's aero model is set to `TAYLOR`. +The samples are the polar, so a stall knee between two of them is represented rather +than smoothed into a slope — which is the whole reason a live source samples instead +of fitting. -`shape` is the [`KulfanParameters`](@ref) the coefficients were generated from, stored -on the panel as `live_shape`. It is the object that was sampled, not a copy or a -re-derivation, and it is written here so a panel's polar and the shape behind it are -set together and cannot drift apart — which is what makes a plot of `live_shape` a -picture of what the solve actually flew rather than of what it should have. +`shape` is the [`KulfanParameters`](@ref) the values were generated from, stored on +the panel as `live_shape` so a panel's polar and the shape behind it are set +together and cannot drift apart. """ -function set_taylor_polar!(panel::Panel, alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs; - window=0.0, shape=nothing) - panel.aero_model = TAYLOR - panel.alpha_ref = Float64(alpha_ref) - panel.alpha_window = Float64(window) +function set_sampled_polar!(panel::Panel, alphas, cl, cd, cm; shape=nothing) + length(alphas) == length(cl) == length(cd) == length(cm) || + throw(ArgumentError("A SAMPLED polar needs one value per angle; got " * + "$(length(alphas)) angles and $(length(cl))/$(length(cd))/" * + "$(length(cm)) values.")) + issorted(alphas) || + throw(ArgumentError("A SAMPLED polar needs ascending angles.")) + panel.aero_model = SAMPLED + panel.alpha_ref = alphas[(length(alphas) + 1) ÷ 2] + panel.alpha_window = maximum(abs, alphas .- panel.alpha_ref) panel.live_shape = shape - for (dst_sym, src) in ((:cl_coeffs, cl_coeffs), (:cd_coeffs, cd_coeffs), - (:cm_coeffs, cm_coeffs)) + for (dst_sym, src) in ((:alpha_knots, alphas), (:cl_coeffs, cl), + (:cd_coeffs, cd), (:cm_coeffs, cm)) dst = getfield(panel, dst_sym) if length(dst) == length(src) dst .= src else - setfield!(panel, dst_sym, Vector{Float64}(src)) + setfield!(panel, dst_sym, collect(Float64, src)) end end return nothing @@ -410,9 +410,8 @@ function calculate_cl(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} cl = 2 * cos(alpha) * sin(alpha)^2 end return R(cl) - elseif panel.aero_model == TAYLOR - return R(taylor_value(panel.cl_coeffs, alpha - panel.alpha_ref, - panel.alpha_window)) + elseif panel.aero_model == SAMPLED + return R(sampled_value(panel.alpha_knots, panel.cl_coeffs, alpha)) elseif panel.aero_model == INVISCID return R(2π * alpha) end @@ -441,12 +440,9 @@ function calculate_cd(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} return R(2 * sin(alpha)^3) end return R(evalpoly(rad2deg(alpha), panel.cd_coeffs)) - elseif panel.aero_model == TAYLOR - # A drag fit has its minimum inside the window, so the slope at the lower - # edge points down and the continuation would carry it through zero into a - # negative drag — an energy source, not an extrapolation. - return R(max(zero(R), taylor_value(panel.cd_coeffs, alpha - panel.alpha_ref, - panel.alpha_window))) + elseif panel.aero_model == SAMPLED + return R(max(zero(R), + sampled_value(panel.alpha_knots, panel.cd_coeffs, alpha))) elseif panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES) cd_interp = panel.cd_interp cd_interp === nothing && @@ -473,9 +469,8 @@ function calculate_cm(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} isnan(alpha) && return R(NaN) if panel.aero_model == POLY return R(evalpoly(rad2deg(alpha), panel.cm_coeffs)) - elseif panel.aero_model == TAYLOR - return R(taylor_value(panel.cm_coeffs, alpha - panel.alpha_ref, - panel.alpha_window)) + elseif panel.aero_model == SAMPLED + return R(sampled_value(panel.alpha_knots, panel.cm_coeffs, alpha)) elseif panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES) cm_interp = panel.cm_interp cm_interp === nothing && diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 3a4c5047..1025ed41 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -1287,20 +1287,6 @@ function calculate_new_aero_data(aero_model, return (alpha_left, delta_left, CL_data, CD_data, CM_data) - elseif isequal(model_type, TAYLOR) - data_left = aero_data[section_index] - data_right = aero_data[section_index + 1] - (data_left isa Tuple{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}} && - data_right isa Tuple{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}}) || - throw(ArgumentError("TAYLOR requires aero_data = " * - "(alpha_ref, cl_coeffs, cd_coeffs, cm_coeffs).")) - return ( - data_left[1] * left_weight + data_right[1] * right_weight, - data_left[2] .* left_weight .+ data_right[2] .* right_weight, - data_left[3] .* left_weight .+ data_right[3] .* right_weight, - data_left[4] .* left_weight .+ data_right[4] .* right_weight, - ) - elseif isequal(model_type, POLY) data_left = aero_data[section_index] data_right = aero_data[section_index + 1] diff --git a/test/airfoil_aero/test_live_polar.jl b/test/airfoil_aero/test_live_polar.jl index 5b443bf6..223b5c7d 100644 --- a/test/airfoil_aero/test_live_polar.jl +++ b/test/airfoil_aero/test_live_polar.jl @@ -4,61 +4,7 @@ using Statistics using VortexStepMethod using VortexStepMethod.AirfoilAero using VortexStepMethod: Panel, calculate_cl, calculate_cd, calculate_cm, - set_taylor_polar! - -@testset "TAYLOR panel polar" begin - panel = Panel{Float64}() - set_taylor_polar!(panel, deg2rad(5.0), [0.6, 5.0, -2.0], [0.02, 0.1, 1.0], - [-0.1, 0.2, 0.0]) - @test panel.aero_model == TAYLOR - alpha = deg2rad(6.0) - d = alpha - deg2rad(5.0) - @test calculate_cl(panel, alpha) ≈ 0.6 + 5.0d - 2.0d^2 - @test calculate_cd(panel, alpha) ≈ 0.02 + 0.1d + 1.0d^2 - @test calculate_cm(panel, alpha) ≈ -0.1 + 0.2d - - # A local expansion has no flap axis: delta must not change the answer. - @test calculate_cl(panel, alpha, deg2rad(9.0)) == calculate_cl(panel, alpha) - - coeffs = panel.cl_coeffs - set_taylor_polar!(panel, 0.0, [1.0, 2.0, 3.0], [1.0, 2.0, 3.0], [1.0, 2.0, 3.0]) - @test panel.cl_coeffs === coeffs # same order refits in place -end - -@testset "TAYLOR window is continued linearly" begin - window = deg2rad(4.0) - panel = Panel{Float64}() - set_taylor_polar!(panel, 0.0, [0.6, 5.0, -20.0], [0.02, 0.0, 0.0], - [-0.1, 0.0, 0.0]; window) - # Inside the window the polynomial is untouched. - @test calculate_cl(panel, 0.5window) ≈ 0.6 + 5.0 * 0.5window - 20.0 * (0.5window)^2 - # At the edge value and slope match, so the continuation is smooth. - edge = calculate_cl(panel, window) - slope = (calculate_cl(panel, window + 1e-7) - edge) / 1e-7 - inner = (edge - calculate_cl(panel, window - 1e-7)) / 1e-7 - @test slope ≈ inner rtol = 1e-4 - # Drag is the exception: its fit has a minimum inside the window, so continuing - # the lower edge's slope would take it through zero. A negative drag coefficient - # is an energy source, and on the SK100 it drove the apparent wind from 12 to - # 21.8 m/s inside one step before the solve died. - drag = Panel{Float64}() - set_taylor_polar!(drag, 0.0, [0.5, 0.0, 0.0], [0.02, 0.4, 6.0], - [-0.05, 0.0, 0.0]; window) - @test calculate_cd(drag, -0.5window) < calculate_cd(drag, 0.0) # still falling - @test calculate_cd(drag, -window) >= 0.0 - @test calculate_cd(drag, -4window) >= 0.0 - @test calculate_cd(drag, -20window) >= 0.0 - # Above the window the continuation is untouched, drag grows. - @test calculate_cd(drag, 4window) > calculate_cd(drag, window) > 0.0 - - # And it stays linear rather than falling off with the parabola's arm, on both - # sides: the unbounded quadratic would be 0.6 + 5·d − 20·d² far out. - @test calculate_cl(panel, 4window) ≈ edge + slope * 3window rtol = 1e-5 - lower_edge = calculate_cl(panel, -window) - lower_slope = 5.0 - 2 * 20.0 * (-window) - @test calculate_cl(panel, -4window) ≈ lower_edge - lower_slope * 3window rtol = 1e-5 - @test calculate_cl(panel, 4window) > 0.6 + 5.0 * 4window - 20.0 * (4window)^2 -end + set_sampled_polar! @testset "Kulfan deformation" begin basis = KulfanBasis() @@ -129,6 +75,44 @@ end @test_throws ArgumentError control_point_deflection(basis, [0.0, 0.0], [1.0, 2.0]) end +@testset "SAMPLED panel polar" begin + panel = Panel{Float64}() + alphas = deg2rad.([-6.0, -3.0, 0.0, 3.0, 6.0]) + set_sampled_polar!(panel, alphas, [0.0, 0.3, 0.6, 0.9, 0.8], + [0.03, 0.02, 0.02, 0.03, 0.06], [-0.1, -0.1, -0.1, -0.1, 0.0]) + @test panel.aero_model == SAMPLED + @test panel.alpha_ref ≈ 0.0 + @test panel.alpha_window ≈ deg2rad(6.0) + + # A sample is the polar at its own angle, and between two it interpolates. + @test calculate_cl(panel, deg2rad(3.0)) ≈ 0.9 + @test calculate_cl(panel, deg2rad(1.5)) ≈ 0.75 + @test calculate_cd(panel, deg2rad(-4.5)) ≈ 0.025 + @test calculate_cm(panel, deg2rad(4.5)) ≈ -0.05 + + # Which is what a fit cannot do: the peak is held rather than averaged away. + @test calculate_cl(panel, deg2rad(6.0)) < calculate_cl(panel, deg2rad(3.0)) + + # Past either end the last sample is held, so a solve that steps out reads a + # value the network gave rather than an extrapolation that runs away. + @test calculate_cl(panel, deg2rad(40.0)) ≈ 0.8 + @test calculate_cl(panel, deg2rad(-40.0)) ≈ 0.0 + @test calculate_cd(panel, deg2rad(-40.0)) ≈ 0.03 + + # A sampled polar has no flap axis: delta must not change the answer. + @test calculate_cl(panel, deg2rad(1.5), deg2rad(9.0)) == + calculate_cl(panel, deg2rad(1.5)) + + knots, coeffs = panel.alpha_knots, panel.cl_coeffs + set_sampled_polar!(panel, alphas .+ deg2rad(2.0), [0.2, 0.5, 0.8, 1.1, 1.0], + [0.03, 0.02, 0.02, 0.03, 0.06], zeros(5)) + @test panel.alpha_knots === knots && panel.cl_coeffs === coeffs + + @test_throws ArgumentError set_sampled_polar!(panel, alphas, [0.0], [0.0], [0.0]) + @test_throws ArgumentError set_sampled_polar!(panel, reverse(alphas), zeros(5), + zeros(5), zeros(5)) +end + @testset "live polars" begin base = KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) n_panels = 3 @@ -136,23 +120,35 @@ end live = LivePolars(fill(base, n_panels)) confidence = refresh_live_polars!(live, panels, deg2rad(6.0), 3e6) @test 0.0 < confidence <= 1.0 - @test all(p -> p.aero_model == TAYLOR, panels) + @test all(p -> p.aero_model == SAMPLED, panels) - # The fit tracks a direct NeuralFoil sweep across its own window. - errors = [abs(calculate_cl(panels[1], deg2rad(6.0 + d)) - - neuralfoil_aero(base, 6.0 + d, 3e6).CL[1]) - for d in range(-4, 4, 17)] - @test maximum(errors) < 0.02 + # Every sample is the network's own answer, not a fit through it. + for (k, offset) in enumerate(live.settings.offsets) + @test calculate_cl(panels[1], deg2rad(6.0) + offset) ≈ + neuralfoil_aero(base, 6.0 + rad2deg(offset), 3e6).CL[1] rtol = 1e-6 + end + # And between them it tracks a direct sweep. Two bounds, because they are two + # different claims: away from the knee the polar is nearly straight and the + # samples pin it down, while at the knee the error is set by how sharp the + # corner is rather than by the spacing — halving the spacing there buys almost + # nothing, and the value stays bounded either way, which is what the solve needs. + sweep(range_deg) = maximum(abs(calculate_cl(panels[1], deg2rad(6.0 + d)) - + neuralfoil_aero(base, 6.0 + d, 3e6).CL[1]) + for d in range_deg) + @test sweep(range(-12, 6, 37)) < 0.02 + @test sweep(range(-12, 12, 49)) < 0.08 @test polar_drift(live, fill(deg2rad(6.0), n_panels)) ≈ 0.0 atol = 1e-12 - @test polar_drift(live, fill(deg2rad(10.0), n_panels)) ≈ 1.0 + @test polar_drift(live, fill(deg2rad(18.0), n_panels)) ≈ 1.0 # Deforming the camber up raises lift at the same angle of attack. flat = calculate_cl(panels[1], deg2rad(6.0)) camber = @. 0.02 * live.basis.x * (1 - live.basis.x) + knots = panels[1].alpha_knots refresh_live_polars!(live, panels, deg2rad(6.0), 3e6; deflection=fill(camber, n_panels)) @test calculate_cl(panels[1], deg2rad(6.0)) > flat + @test panels[1].alpha_knots === knots # a refresh writes in place # The panel keeps the object that was sampled, not a copy of it: a plot of # live_shape has to be a picture of what the solve flew. @@ -165,21 +161,101 @@ end @test_throws ArgumentError refresh_live_polars!(live, panels[1:2], 0.0, 3e6) @test_throws ArgumentError LivePolars(fill(base, 2); - settings=LivePolarSettings(; order=4, n_samples=3)) + settings=LivePolarSettings(; offsets=deg2rad.([1.0, 2.0, 3.0]))) + @test_throws ArgumentError LivePolars(fill(base, 2); + settings=LivePolarSettings(; offsets=deg2rad.([4.0, -4.0]))) end -@testset "TAYLOR spanwise blend" begin - left = (deg2rad(4.0), [0.5, 5.0, 0.0], [0.02, 0.0, 0.0], [-0.1, 0.0, 0.0]) - right = (deg2rad(8.0), [0.9, 5.0, 0.0], [0.06, 0.0, 0.0], [-0.3, 0.0, 0.0]) - blended = VortexStepMethod.calculate_new_aero_data( - (TAYLOR, TAYLOR), (left, right), 1, 0.25, 0.75) - @test blended[1] ≈ deg2rad(7.0) - @test blended[2] ≈ [0.8, 5.0, 0.0] +@testset "live surface pressure" begin + base = KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) + n_panels = 2 + panels = [Panel{Float64}() for _ in 1:n_panels] + settings = LivePolarSettings(; model_size="large") + live = LivePolars(fill(base, n_panels); settings) + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6) - section_left = Section([0.0, 1.0, 0.0], [1.0, 1.0, 0.0], TAYLOR, left) - section_right = Section([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], TAYLOR, right) - panel = Panel{Float64}() - VortexStepMethod.init_aero!(panel, section_left, section_right) - @test panel.alpha_ref ≈ deg2rad(6.0) - @test panel.cl_coeffs ≈ [0.7, 5.0, 0.0] + # The contour the traction pattern is spread over: Selig order, nose at argmin. + contour_x, contour_y = kulfan_to_coordinates(base) + contour_x = collect(float.(contour_x)) + leading_edge = argmin(contour_x) + cp = [zeros(length(contour_x)) for _ in 1:n_panels] + same = cp[1] + refresh_live_pressure!(cp, live, fill(contour_x, n_panels), + fill(leading_edge, n_panels), deg2rad(6.0), 3e6) + @test cp[1] === same # written in place + @test all(isfinite, cp[1]) + @test cp[1] ≈ cp[2] # same shape, same angle + + # This is the same pressure the offline tables were generated with, so the live + # path and the table generator must not be two different answers for one shape. + solver = NeuralFoilSolver(; model_size="large", n_crit=settings.n_crit) + section = analyze_section(solver, DeformedSection(base, contour_x, + collect(float.(contour_y))), deg2rad(6.0), 3e6) + @test cp[1] ≈ section.cp rtol = 1e-8 + + # Suction over the upper surface, pressure under the lower one. + @test minimum(cp[1][1:leading_edge]) < -0.5 + loading(p) = sum(view(p, (leading_edge + 1):length(p))) / + (length(p) - leading_edge) - sum(view(p, 1:leading_edge)) / + leading_edge + @test loading(cp[1]) > 0.0 + + # Bending the camber up loads the section harder, which is the whole claim: the + # pattern moves with the deformation. Read as net loading, not as peak suction — + # a cambered nose meets the flow better, so the leading-edge peak actually + # weakens (-2.64 to -2.57) while the section as a whole carries more. + camber = @. 0.02 * live.basis.x * (1 - live.basis.x) + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6; + deflection=fill(camber, n_panels)) + cambered = deepcopy(cp) + refresh_live_pressure!(cambered, live, fill(contour_x, n_panels), + fill(leading_edge, n_panels), deg2rad(6.0), 3e6) + @test loading(cambered[1]) > loading(cp[1]) + + @test_throws ArgumentError refresh_live_pressure!(cp[1:1], live, + fill(contour_x, n_panels), fill(leading_edge, n_panels), 0.0, 3e6) +end + +@testset "the contour follows the deformation" begin + base = KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) + panels = [Panel{Float64}() for _ in 1:2] + live = LivePolars(fill(base, 2)) + shape = [contour_shape_matrix(live.basis.x, 8) for _ in 1:2] + offset = [zeros(length(live.basis.x)) for _ in 1:2] + same = offset[1] + + # An undeformed panel's contour is the reference one, exactly. + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6) + live_shape_offset!(offset, live, shape) + @test offset[1] === same + @test all(iszero, offset[1]) + + # A camber increment comes back as itself: this is the same offset the airfoil + # the network saw was built from, so the contour and the pressures agree. + camber = @. 0.02 * live.basis.x * (1 - live.basis.x) + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6; + deflection=fill(camber, 2)) + live_shape_offset!(offset, live, shape) + @test maximum(abs, offset[1] .- camber) < 0.02 * maximum(abs, camber) + @test offset[1] ≈ offset[2] + + # A chord rotation is not a shape change, so it moves no contour node. + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6; + deflection=fill(0.05 .* live.basis.x, 2)) + live_shape_offset!(offset, live, shape) + @test maximum(abs, offset[1]) < 1e-12 + + @test_throws ArgumentError live_shape_offset!(offset[1:1], live, shape) +end + +@testset "live skin friction" begin + contour_x = [1.0, 0.5, 0.0, 0.5, 1.0] + cf = [zeros(5), zeros(5)] + same = cf[1] + live_surface_friction!(cf, [contour_x, contour_x], [3e6, 1e6]) + @test cf[1] === same + @test cf[1] ≈ [AirfoilAero.flat_plate_cf(x, 3e6) for x in contour_x] + # Lower Reynolds is more friction, and it is highest at the nose. + @test all(cf[2] .> cf[1]) + @test argmax(cf[1]) == 3 end diff --git a/test/plotting/test_plotting.jl b/test/plotting/test_plotting.jl index e3c5dbc8..afe3c9ba 100644 --- a/test/plotting/test_plotting.jl +++ b/test/plotting/test_plotting.jl @@ -528,8 +528,9 @@ end cambered = VortexStepMethod.AirfoilAero.deform_kulfan(basis, base, @. 0.08 * basis.x * (1 - basis.x)) for panel in body_aero.panels - VortexStepMethod.set_taylor_polar!(panel, 0.0, [0.5, 5.0, 0.0], - [0.02, 0.0, 0.0], [-0.05, 0.0, 0.0]; window=deg2rad(4), shape=cambered) + VortexStepMethod.set_sampled_polar!(panel, deg2rad.([-4.0, 0.0, 4.0]), + [0.3, 0.5, 0.7], [0.02, 0.02, 0.03], [-0.05, -0.05, -0.05]; + shape=cambered) end @test body_aero.panels[1].live_shape === cambered v_live, f_live, ribs_live = airfoil_skin_geometry(body_aero)