From 1e8f13367cb7d01aa8fce18ce2825499dd36437f Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 31 Aug 2026 20:38:07 +0200 Subject: [PATCH 1/2] A live polar is a polar, not a model of its own SAMPLED was a second linear interpolation over a second set of knots, doing what POLAR_VECTORS already does with Interpolations, and diverging from it: drag extrapolated linearly under one and was held flat under the other, and delta needed a rule of its own. What a live source actually needs is not a model, it is the ability to rewrite a table it already has. set_polar! rewrites a POLAR_VECTORS panel's knots and values in its own storage and rebuilds the interpolations, since Interpolations copies the values it is handed. The panel records the range its table covers, and window_alpha holds alpha inside it, so a table sampled over a window around one angle stops at its ends instead of extrapolating out of them. That is a property of a bounded table rather than of a live one, so the stall scan can key off it and full-range polars, whose window is zero, are untouched. init_aero! now clears the window and the live shape, so a mesh rebuild that re-seeds a panel from its section cannot leave a stale window clamping the section's own table. SAMPLED never shipped, so nothing here is breaking. Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 60 +++++++-------- docs/src/private_functions.md | 5 +- src/VortexStepMethod.jl | 18 ++--- src/airfoil_aero/AirfoilAero.jl | 2 +- src/airfoil_aero/live_polar.jl | 16 ++-- src/body_aerodynamics.jl | 7 +- src/panel.jl | 105 ++++++++++++++------------- test/airfoil_aero/test_live_polar.jl | 46 ++++++++---- test/plotting/test_plotting.jl | 9 ++- 9 files changed, 145 insertions(+), 123 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 48c48615..0a4d3a12 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,7 +10,7 @@ nothing over the first and last `1/2N` of chord; the old code extrapolated `Cp` off the end of each surface independently. Through a stagnation point `ue` is linear in arc length while `Cp` is quadratic, so that extrapolated a - parabola through its own turning point: on an SK100 section it put + parabola through its own turning point: on an inflated section it put `Cp = -1.60` at the leading edge, where it has to approach `+1`, and the two surfaces disagreed at the trailing edge in violation of the Kutta condition. Signing the lower surface negative makes both surfaces one continuous curve @@ -22,29 +22,31 @@ ### Added -- BREAKING: the `TAYLOR` aero model is gone, with `taylor_value` and - `set_taylor_polar!`. An order-2 polynomial cannot hold a stall knee: fitted - over a window wider than the knee it averages across it, narrower and it never - reaches it, and past the window it continued on whatever edge slope it ended - with — which on an SK100 tip panel was negative where the local slope was - `+0.05`/deg, and took the circulation solve to `NaN`. `SAMPLED` replaces it. -- `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. +- A `POLAR_VECTORS` panel's table can be rewritten at run time by `set_polar!`, + which is what a live polar source does every solve. The knots and values reuse + the panel's own storage, and the panel records the range the table covers, so + a table generated over a window around one angle of attack is held at its end + values rather than extrapolated past them, and the stall-angle scan skips a + table that stops short of it. A rewritten table is an ordinary polar in every + other way: same evaluation, same `delta` handling, no separate aero model. - 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 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. + onto that basis. `chord_residual` first takes the straight line through the + deflection's own endpoints off it, which the basis cannot express and answers + with runaway weights; `chord_line` returns that line, since its rotation is + one the caller may still owe its angle of attack. `panel_kulfan_parameters` + fits each panel's contour once at build time, shrink-wrapped first because a + raw slice fits to weights that oscillate far more than any deformation will + and leave neighbouring panels disagreeing enough to cost the solve its + convergence. `LivePolars` and `refresh_live_polars!` then evaluate NeuralFoil + on a grid of angles per panel in one batched forward pass and write those + values in as each panel's polar table, 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 @@ -61,15 +63,15 @@ table generator reconstruct pressure the same way. - `compare_live_polar` solves one panel's current deformed shape in XFoil and puts the answer next to the live polar the panel is flying, at the same angle - and Reynolds and under `live_xfoil_solver`'s matched transition settings. XFoil - is marched out from zero rather than jumped to the angle — the viscous march - refuses a blunt inflated section otherwise, and on the SK100's deformed - mid-span it refused every one — and the comparison is made at the nearest angle - it reached, reported as `alpha` beside the `requested` one. A - `LivePolars` now keeps each panel's analysis confidence, so the panel worth - asking about is the one to hand. NeuralFoil's confidence scores the input - shape rather than the answer, so a deformed section far from its training set - reads low whether or not the polar is wrong; this is what tells the two apart. + and Reynolds and under `live_xfoil_solver`'s matched transition settings. + XFoil is marched out from zero rather than jumped to the angle — the viscous + march refuses a blunt inflated section otherwise, and on a deformed mid-span + panel it refused every one — and the comparison is made at the nearest angle + it reached, reported as `alpha` beside the `requested` one. A `LivePolars` now + keeps each panel's analysis confidence, so the panel worth asking about is the + one to hand. NeuralFoil's confidence scores the input shape rather than the + answer, so a deformed section far from its training set reads low whether or + not the polar is wrong; this is what tells the two apart. - `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 diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 0f579ab2..1bc9c242 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -98,8 +98,9 @@ update_non_deformed_sections! ### Aerodynamic data and Cp ```@docs calculate_new_aero_data -set_sampled_polar! -sampled_value +set_polar! +rebuild_polar +window_alpha assemble_polar_matrix load_matrix_polar_data read_aero_matrix diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index badf43b3..375d664f 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, SAMPLED +export AeroModel, INVISCID, POLY, LEI_AIRFOIL_BREUKELS, POLAR_MATRICES, POLAR_VECTORS export KulfanParameters export BILLOWING, COSINE, LINEAR, PanelDistribution, SPLIT_PROVIDED, UNCHANGED export ELLIPTIC, InitialGammaDistribution, ZEROS @@ -238,21 +238,19 @@ Enumeration of the implemented wing types. @enum WingType RECTANGULAR CURVED ELLIPTICAL """ - AeroModel `POLY` `POLAR_VECTORS` `POLAR_MATRICES` `INVISCID` `SAMPLED` + AeroModel `POLY` `POLAR_VECTORS` `POLAR_MATRICES` `INVISCID` Enumeration of the implemented aerodynamic models. See also: [`AeroData`](@ref) # Elements - `POLY`: α-polynomial coefficients for cl/cd/cm (e.g. Breukels LEI coeffs, generated by the `AirfoilAero` package). Core only evaluates the polynomial. -- `POLAR_VECTORS`: Polar vectors as function of alpha (lookup tables with interpolation) +- `POLAR_VECTORS`: Polar vectors as function of alpha (lookup tables with interpolation). + A panel's table may be rewritten at run time by [`set_polar!`](@ref + VortexStepMethod.set_polar!), which is how a live polar source regenerates it from the + panel's deformed shape each solve. - `POLAR_MATRICES`: Polar matrices as function of alpha and delta (lookup tables with interpolation) - INVISCID -- `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`. @@ -263,7 +261,6 @@ where `alpha` is the angle of attack, `delta` is trailing edge angle. POLAR_VECTORS POLAR_MATRICES INVISCID - SAMPLED end """ @@ -346,9 +343,6 @@ Union of different definitions of the aerodynamic properties of a wing section. - (`alpha_range`, `cl_vector`, `cd_vector`, `cm_vector`) for `POLAR_VECTORS` - (`alpha_range`, `delta_range`, `cl_matrix`, `cd_matrix`, `cm_matrix`) for `POLAR_MATRICES` -`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 the curvature of its airfoil shape. The camber is typically measured as the maximum distance diff --git a/src/airfoil_aero/AirfoilAero.jl b/src/airfoil_aero/AirfoilAero.jl index dbc14589..140da02a 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_sampled_polar!, + write_node_rows, section_surface, set_polar!, KulfanParameters, calculate_cl, calculate_cd, calculate_cm diff --git a/src/airfoil_aero/live_polar.jl b/src/airfoil_aero/live_polar.jl index 733af75d..215eccbc 100644 --- a/src/airfoil_aero/live_polar.jl +++ b/src/airfoil_aero/live_polar.jl @@ -4,7 +4,7 @@ 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. +panel's polar table directly — there is no fit in between. 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 @@ -156,14 +156,15 @@ end refresh_live_polars!(live, panels, alpha_ref, reynolds; deflection=nothing) -> Float64 -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: +Regenerate every panel's polar table from its current shape and write it in +(see [`set_polar!`](@ref VortexStepMethod.set_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. +the forward pass and the three interpolation rebuilds `Interpolations` needs to take the +new values, which it copies rather than references. Each panel keeps the deformed shape it was evaluated at as its `live_shape`, so a plot draws the airfoil the network actually saw. @@ -204,9 +205,8 @@ function refresh_live_polars!(live::LivePolars, panels, alpha_ref, reynolds; samples = ((i - 1) * n_samples + 1):(i * n_samples) live.knots .= alpha_vec[i] .+ offsets live.confidence[i] = minimum(view(confidence, samples)) - set_sampled_polar!(panels[i], live.knots, view(cl, samples), - view(cd, samples), view(cm, samples); - shape=live.deformed[i]) + set_polar!(panels[i], live.knots, view(cl, samples), view(cd, samples), + view(cm, samples); shape=live.deformed[i]) end return minimum(confidence) end @@ -239,7 +239,7 @@ pair stays like for like even when the march stops short. Run it when `confidence` is low. A confidence is the network's opinion of its own inputs — a shape far from what it was trained on scores badly whether or not the answer is wrong — so it says to go and check, not what the check will find. `cl` here is read -off the panel's `SAMPLED` polar rather than from a fresh network call, so what is +off the panel's polar table rather than from a fresh network call, so what is compared is what the solver actually flew. A reference that will not converge anywhere on the ramp comes back `NaN` rather than diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 16d8b7c0..8ecd60ff 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -224,9 +224,10 @@ function calculate_stall_angle_list!(stall_angles::AbstractVector, # Default stall angle if none found panel_stall = stall_angle_if_none_detected - # 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 + # A table that stops short of the scan says nothing about a stall in it: + # past its end it is held flat, which is not a peak. + if panel.alpha_window > 0 && + panel.alpha_ref + panel.alpha_window < deg2rad(begin_aoa) stall_angles[idx] = panel_stall continue end diff --git a/src/panel.jl b/src/panel.jl index b47014fa..63b6f127 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -36,9 +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: 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 +- `alpha_ref`::Float64=0: centre angle [rad] of the range the polar table covers +- `alpha_window`::Float64=0: half width [rad] the table reaches; 0 = unbounded +- `alpha_knots`::Vector{Float64}=Float64[]: ascending angles [rad] the panel's own table 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} @@ -214,6 +214,8 @@ end function init_aero!(panel::Panel, section_1::Section, section_2::Section; remove_nan = true) panel.aero_model = section_1.aero_model + panel.alpha_ref, panel.alpha_window = 0.0, 0.0 + panel.live_shape = nothing section_1.aero_model == section_2.aero_model || throw(ArgumentError("Both sections must have the same aero model, not " * "$(section_1.aero_model) and $(section_2.aero_model)")) @@ -233,61 +235,54 @@ function init_aero!(panel::Panel, section_1::Section, section_2::Section; end """ - sampled_value(knots, values, alpha) + window_alpha(panel, alpha) -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. - -`knots` is ascending; a single knot is a constant polar. +`alpha` [rad] clamped into the range the panel's polar was built over. A table generated +over a window around one angle of attack says nothing past its ends, so it is held at its +end values rather than extrapolated out of them; `alpha_window` of `0` means unbounded and +leaves `alpha` alone, which is what a full-range polar wants. """ -@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 +@inline function window_alpha(panel::Panel, alpha) + panel.alpha_window > 0 || return alpha + return clamp(alpha, panel.alpha_ref - panel.alpha_window, + panel.alpha_ref + panel.alpha_window) end +set_polar!(panel::Panel{<:Any, Nothing}, alphas, cl, cd, cm; shape=nothing) = + throw(ArgumentError("set_polar! needs a panel that carries interpolations. This one " * + "was built without them (INVISCID or POLY sections); build the wing from " * + "POLAR_VECTORS sections so its panels have the interpolation types.")) + """ - set_sampled_polar!(panel, alphas, cl, cd, cm; shape=nothing) + set_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`. +Rewrite a panel's `POLAR_VECTORS` table with values at `alphas` [rad], ascending, and +rebuild its interpolations. The panel's `alpha_ref` and `alpha_window` are taken from the +angles, so a table covering only a window around one angle is held at its ends rather than +extrapolated past them (see [`window_alpha`](@ref)). 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. - -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. +`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 growing the +panel. The interpolations themselves are rebuilt, since `Interpolations` copies the values +it is handed; they keep the extrapolation the panel was built with. -`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. +`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_sampled_polar!(panel::Panel, alphas, cl, cd, cm; shape=nothing) +function set_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 " * + throw(ArgumentError("A polar table needs one value per angle; got " * "$(length(alphas)) angles and $(length(cl))/$(length(cd))/" * "$(length(cm)) values.")) + length(alphas) >= 2 || + throw(ArgumentError("A polar table needs at least two angles.")) 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) + throw(ArgumentError("A polar table needs ascending angles.")) + panel.aero_model = POLAR_VECTORS + panel.alpha_ref = (alphas[1] + alphas[end]) / 2 + panel.alpha_window = (alphas[end] - alphas[1]) / 2 panel.live_shape = shape for (dst_sym, src) in ((:alpha_knots, alphas), (:cl_coeffs, cl), (:cd_coeffs, cd), (:cm_coeffs, cm)) @@ -298,9 +293,21 @@ function set_sampled_polar!(panel::Panel, alphas, cl, cd, cm; shape=nothing) setfield!(panel, dst_sym, collect(Float64, src)) end end + panel.cl_interp = rebuild_polar(panel.cl_interp, panel.alpha_knots, panel.cl_coeffs) + panel.cd_interp = rebuild_polar(panel.cd_interp, panel.alpha_knots, panel.cd_coeffs) + panel.cm_interp = rebuild_polar(panel.cm_interp, panel.alpha_knots, panel.cm_coeffs) return nothing end +""" + rebuild_polar(old, knots, values) -> Extrapolation + +A 1D linear interpolation over `knots`/`values` carrying `old`'s extrapolation, so the +rebuilt object has the type the panel's field was parameterised with. +""" +rebuild_polar(old, knots, values) = + linear_interpolation(knots, values; extrapolation_bc=old.et) + """ 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...) @@ -404,14 +411,13 @@ calculate_cl(panel::Panel, alpha) = calculate_cl(panel, alpha, panel.delta) function calculate_cl(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} R = promote_type(Tp, Ta, Td) isnan(alpha) && return R(NaN) + alpha = window_alpha(panel, alpha) if panel.aero_model == POLY cl = evalpoly(rad2deg(alpha), panel.cl_coeffs) if abs(alpha) > (π/9) cl = 2 * cos(alpha) * sin(alpha)^2 end return R(cl) - 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 @@ -435,14 +441,12 @@ calculate_cd(panel::Panel, alpha) = calculate_cd(panel, alpha, panel.delta) function calculate_cd(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} R = promote_type(Tp, Ta, Td) isnan(alpha) && return R(NaN) + alpha = window_alpha(panel, alpha) if panel.aero_model == POLY if abs(alpha) > (π/9) # Outside ±20 degrees return R(2 * sin(alpha)^3) end return R(evalpoly(rad2deg(alpha), panel.cd_coeffs)) - 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 && @@ -467,10 +471,9 @@ calculate_cm(panel::Panel, alpha) = calculate_cm(panel, alpha, panel.delta) function calculate_cm(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} R = promote_type(Tp, Ta, Td) isnan(alpha) && return R(NaN) + alpha = window_alpha(panel, alpha) if panel.aero_model == POLY return R(evalpoly(rad2deg(alpha), panel.cm_coeffs)) - 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/test/airfoil_aero/test_live_polar.jl b/test/airfoil_aero/test_live_polar.jl index 70bd8878..dafaedc3 100644 --- a/test/airfoil_aero/test_live_polar.jl +++ b/test/airfoil_aero/test_live_polar.jl @@ -3,8 +3,21 @@ using LinearAlgebra using Statistics using VortexStepMethod using VortexStepMethod.AirfoilAero -using VortexStepMethod: Panel, calculate_cl, calculate_cd, calculate_cm, - set_sampled_polar! +using VortexStepMethod: Panel, Section, calculate_cl, calculate_cd, calculate_cm, + set_polar!, panel_interp_types, init_aero! + +"""Panels carrying the interpolations a `POLAR_VECTORS` table is rewritten into.""" +function polar_panels(n_panels; n_samples=5) + alphas = collect(range(deg2rad(-6.0), deg2rad(6.0), n_samples)) + data = (alphas, zeros(n_samples), fill(0.01, n_samples), zeros(n_samples)) + section = Section([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], POLAR_VECTORS, data) + CL, CD, CM, CP = panel_interp_types(section, true) + return map(1:n_panels) do _ + panel = Panel{Float64, CL, CD, CM, CP}() + init_aero!(panel, section, section) + panel + end +end @testset "Kulfan deformation" begin basis = KulfanBasis() @@ -75,12 +88,12 @@ 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}() +@testset "rewritten panel polar" begin + panel = only(polar_panels(1)) 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], + set_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.aero_model == POLAR_VECTORS @test panel.alpha_ref ≈ 0.0 @test panel.alpha_window ≈ deg2rad(6.0) @@ -104,23 +117,28 @@ end 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], + set_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)) + @test_throws ArgumentError set_polar!(panel, alphas, [0.0], [0.0], [0.0]) + @test_throws ArgumentError set_polar!(panel, reverse(alphas), zeros(5), + zeros(5), zeros(5)) + @test_throws ArgumentError set_polar!(panel, alphas[1:1], [0.5], [0.02], [-0.1]) + + # A panel built without interpolations says so rather than failing on a field type. + @test_throws ArgumentError set_polar!(Panel{Float64}(), 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 - panels = [Panel{Float64}() for _ in 1:n_panels] + panels = polar_panels(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 == SAMPLED, panels) + @test all(p -> p.aero_model == POLAR_VECTORS, panels) # Every sample is the network's own answer, not a fit through it. for (k, offset) in enumerate(live.settings.offsets) @@ -179,7 +197,7 @@ end @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] + panels = polar_panels(n_panels) settings = LivePolarSettings(; model_size="large") live = LivePolars(fill(base, n_panels); settings) refresh_live_polars!(live, panels, deg2rad(6.0), 3e6) @@ -232,7 +250,7 @@ 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] + panels = polar_panels(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] diff --git a/test/plotting/test_plotting.jl b/test/plotting/test_plotting.jl index afe3c9ba..b6a973ae 100644 --- a/test/plotting/test_plotting.jl +++ b/test/plotting/test_plotting.jl @@ -481,9 +481,12 @@ function create_body_aero_with_skin(; n_panels=4) cf = fill(0.003, n_node, length(alpha_range), 1) section_aero = SectionAero(alpha_range, [0.0], x, y, cp, cf) + # POLAR_VECTORS, so the panels carry the interpolations a live polar rewrites. + polar = (alpha_range, [0.0, 0.5, 1.0], fill(0.02, 3), fill(-0.05, 3)) wing = Wing(n_panels, spanwise_distribution=LINEAR) - add_section!(wing, [0.0, 2.0, 0.0], [1.0, 2.0, 0.0], INVISCID, nothing, section_aero) - add_section!(wing, [0.0, -2.0, 0.0], [1.0, -2.0, 0.0], INVISCID, nothing, section_aero) + add_section!(wing, [0.0, 2.0, 0.0], [1.0, 2.0, 0.0], POLAR_VECTORS, polar, section_aero) + add_section!(wing, [0.0, -2.0, 0.0], [1.0, -2.0, 0.0], POLAR_VECTORS, polar, + section_aero) refine!(wing) body_aero = BodyAerodynamics([wing]) set_va!(body_aero, [20.0, 0.0, 1.0]) @@ -528,7 +531,7 @@ end cambered = VortexStepMethod.AirfoilAero.deform_kulfan(basis, base, @. 0.08 * basis.x * (1 - basis.x)) for panel in body_aero.panels - VortexStepMethod.set_sampled_polar!(panel, deg2rad.([-4.0, 0.0, 4.0]), + VortexStepMethod.set_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 From 69f74089fe139fa0c3d5b88dbd200bcc003ab55c Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 31 Aug 2026 21:29:26 +0200 Subject: [PATCH 2/2] Let a short polar table scan for its knot instead of bisecting MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Interpolations finds a gridded knot with searchsortedfirst, which dispatches on the knot vector, so the search is ours to choose without changing the library. A binary search over nine angles costs more than the whole lookup around it — its branches are unpredictable where a scan over contiguous memory is not — and the live polars a solve rewrites every step are exactly that short. ScanKnots wraps the angles and scans them, and polar_knots hands it to a table under SCAN_KNOT_MAX and the plain vector to anything longer, where the scan would lose badly: 71 ns against 13 at 361 knots. On the 44 panel wing the solve goes from 215 to 176 us, against 206 for the hand-rolled polar this replaces. Co-Authored-By: Claude Opus 5 --- docs/src/private_functions.md | 2 ++ docs/src/private_types.md | 2 ++ src/panel.jl | 56 +++++++++++++++++++++++++++++++---- 3 files changed, 54 insertions(+), 6 deletions(-) diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 1bc9c242..4979a91f 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -100,6 +100,8 @@ update_non_deformed_sections! calculate_new_aero_data set_polar! rebuild_polar +same_knots +polar_knots window_alpha assemble_polar_matrix load_matrix_polar_data diff --git a/docs/src/private_types.md b/docs/src/private_types.md index ff63c1fc..ff99a2c5 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 +ScanKnots KulfanParameters PanelProperties Filament @@ -17,6 +18,7 @@ SemiInfiniteFilament ### Aerodynamic model constants ```@docs LEI_AIRFOIL_BREUKELS +SCAN_KNOT_MAX ``` ### Airfoil aerodynamics (AirfoilAero) diff --git a/src/panel.jl b/src/panel.jl index 63b6f127..d31d0535 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -132,6 +132,42 @@ function init_pos!( return nothing end +""" + ScanKnots(data) + +Polar angles that a lookup scans rather than bisects. `Interpolations` finds a gridded +knot with `searchsortedfirst`, which dispatches on the knot vector, so a short table can +choose the search that suits it: below `SCAN_KNOT_MAX` a linear scan over contiguous +memory beats a binary search's mispredicted branches, and above it loses badly. Wraps +its vector rather than copying it, so the angles stay writable in place. +""" +struct ScanKnots <: AbstractVector{Float64} + data::Vector{Float64} +end + +"Longest table a linear knot scan still beats a binary search on." +const SCAN_KNOT_MAX = 24 + +Base.size(knots::ScanKnots) = size(knots.data) +Base.@propagate_inbounds Base.getindex(knots::ScanKnots, i::Int) = knots.data[i] +Base.IndexStyle(::Type{ScanKnots}) = IndexLinear() + +@inline function Base.searchsortedfirst(knots::ScanKnots, x, ::Base.Order.ForwardOrdering) + @inbounds for i in eachindex(knots.data) + knots.data[i] >= x && return i + end + return length(knots.data) + 1 +end + +""" + polar_knots(alphas) -> AbstractVector + +The angles of a 1D polar in the container whose knot search fits its length, see +[`ScanKnots`](@ref). +""" +polar_knots(alphas::Vector{Float64}) = + length(alphas) <= SCAN_KNOT_MAX ? ScanKnots(alphas) : alphas + """ build_interps(section_1, section_2, remove_nan) -> (cl, cd, cm, section_aero) @@ -162,9 +198,10 @@ function build_interps(section_1::Section, section_2::Section, remove_nan) cl = Vector{Float64}((aero_1[2] .+ aero_2[2]) ./ 2) cd = Vector{Float64}((aero_1[3] .+ aero_2[3]) ./ 2) cm = Vector{Float64}((aero_1[4] .+ aero_2[4]) ./ 2) - cl_i = linear_interpolation(alphas, cl; extrapolation_bc=extrap_flat) - cd_i = linear_interpolation(alphas, cd; extrapolation_bc=extrap_line) - cm_i = linear_interpolation(alphas, cm; extrapolation_bc=extrap_flat) + knots = polar_knots(alphas) + cl_i = linear_interpolation(knots, cl; extrapolation_bc=extrap_flat) + cd_i = linear_interpolation(knots, cd; extrapolation_bc=extrap_line) + cm_i = linear_interpolation(knots, cm; extrapolation_bc=extrap_flat) else (aero_1 isa Tuple{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}} && aero_2 isa Tuple{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}) || @@ -302,11 +339,18 @@ end """ rebuild_polar(old, knots, values) -> Extrapolation -A 1D linear interpolation over `knots`/`values` carrying `old`'s extrapolation, so the -rebuilt object has the type the panel's field was parameterised with. +A 1D linear interpolation over `knots`/`values` carrying `old`'s extrapolation and knot +container, so the rebuilt object has the type the panel's field was parameterised with. +Only the values need rebuilding — `Interpolations` holds the knots by reference — but +they are what it copies, so the object is rebuilt until it can take them in place. """ rebuild_polar(old, knots, values) = - linear_interpolation(knots, values; extrapolation_bc=old.et) + linear_interpolation(same_knots(old.itp.knots[1], knots), values; + extrapolation_bc=old.et) + +"The angles in the container the panel's interpolations were parameterised with." +same_knots(::ScanKnots, knots) = ScanKnots(knots) +same_knots(::AbstractVector, knots) = knots """ reinit!(panel, section_1, section_2, aero_center, control_point, bound_point_1,