From dad87ab220aae2476f2e32b2c7af3094c812f1e1 Mon Sep 17 00:00:00 2001 From: Bart Date: Tue, 1 Sep 2026 09:30:42 +0200 Subject: [PATCH 1/5] Let a polar refresh write into the panel instead of rebuilding it interpolate! hands back a gridded interpolation that reads the arrays it was given rather than copies of them, knots and values alike, which is what the rebuild in set_polar! was there to work around. Once a panel's interpolations read its own alpha_knots and coefficient vectors, writing those vectors is the whole update: a refresh of an unchanged shape allocates nothing and leaves the objects in place. A table that changes shape is still built once, and the two dimensional form fills every delta column of the matrix it already owns. The assignment back into the panel is guarded, since storing a returned immutable into a mutable field boxes it even when nothing changed. Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 12 +++++++ docs/src/private_functions.md | 3 +- src/panel.jl | 50 +++++++++++++++++++--------- test/airfoil_aero/test_live_polar.jl | 12 +++++++ 4 files changed, 61 insertions(+), 16 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0596fbfd..a072d783 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,17 @@ # Changelog +## Unreleased + +### Changed + +- A live polar refresh writes into the panel's own storage instead of rebuilding + its interpolations. `Interpolations` hands back a gridded interpolation that + reads the very arrays it was given when it is built with `interpolate!`, so + once a panel's table reads its own knots and values, writing them is the whole + update: `set_polar!` on an unchanged table shape allocates nothing and leaves + the interpolation objects in place. A table that changes shape is rebuilt + once. + ## VortexStepMethod v4.3.0 2026-08-31 ### Fixed diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index a9ef5612..fa41f393 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -99,8 +99,9 @@ update_non_deformed_sections! ```@docs calculate_new_aero_data set_polar! -rebuild_polar +refresh_polar! polar_model +reads_from same_knots polar_knots window_alpha diff --git a/src/panel.jl b/src/panel.jl index 3565ce07..2c4de856 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -333,31 +333,51 @@ function set_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) + cl_interp = refresh_polar!(panel.cl_interp, panel.alpha_knots, panel.cl_coeffs) + cd_interp = refresh_polar!(panel.cd_interp, panel.alpha_knots, panel.cd_coeffs) + cm_interp = refresh_polar!(panel.cm_interp, panel.alpha_knots, panel.cm_coeffs) + cl_interp === panel.cl_interp || (panel.cl_interp = cl_interp) + cd_interp === panel.cd_interp || (panel.cd_interp = cd_interp) + cm_interp === panel.cm_interp || (panel.cm_interp = cm_interp) return nothing end """ - rebuild_polar(old, knots, values) -> Extrapolation + refresh_polar!(old, knots, values) -> Extrapolation -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. +An interpolation reading `knots` and `values` themselves, carrying `old`'s extrapolation +and knot container so it keeps the type the panel's field was parameterised with. +`interpolate!` takes both arrays by reference, so once an interpolation reads a panel's +own storage, writing that storage is the whole update and there is nothing to rebuild. +A table that changes shape is rebuilt once, and reads in place from then on. The two +dimensional form spreads the angles over every `delta` the panel spans. """ -rebuild_polar(old::Interpolations.Extrapolation{<:Any, 1}, knots, values) = - linear_interpolation(same_knots(old.itp.knots[1], knots), values; - extrapolation_bc=old.et) +function refresh_polar!(old::E, knots, values + ) where {E <: Interpolations.Extrapolation{<:Any, 1}} + reads_from(old.itp.knots[1], knots) && old.itp.coefs === values && return old + return extrapolate(interpolate!((same_knots(old.itp.knots[1], knots),), values, + Gridded(Linear())), old.et)::E +end -function rebuild_polar(old::Interpolations.Extrapolation{<:Any, 2}, knots, values) +function refresh_polar!(old::E, knots, values + ) where {E <: Interpolations.Extrapolation{<:Any, 2}} deltas = old.itp.knots[2] - return linear_interpolation((same_knots(old.itp.knots[1], knots), deltas), - repeat(values, 1, length(deltas)); - extrapolation_bc=old.et) + coefs = old.itp.coefs + if reads_from(old.itp.knots[1], knots) && size(coefs) == (length(knots), length(deltas)) + for column in axes(coefs, 2) + @views coefs[:, column] .= values + end + return old + end + return extrapolate(interpolate!((same_knots(old.itp.knots[1], knots), deltas), + repeat(values, 1, length(deltas)), + Gridded(Linear())), old.et)::E end +"Whether an interpolation's knot container is a view of `knots` rather than a copy." +reads_from(container::ScanKnots, knots) = container.data === knots +reads_from(container::AbstractVector, knots) = container === knots + "The angles in the container the panel's interpolations were parameterised with." same_knots(::ScanKnots, knots) = ScanKnots(knots) same_knots(::AbstractVector, knots) = knots diff --git a/test/airfoil_aero/test_live_polar.jl b/test/airfoil_aero/test_live_polar.jl index 58afab15..61c9e869 100644 --- a/test/airfoil_aero/test_live_polar.jl +++ b/test/airfoil_aero/test_live_polar.jl @@ -121,6 +121,18 @@ end [0.03, 0.02, 0.02, 0.03, 0.06], zeros(5)) @test panel.alpha_knots === knots && panel.cl_coeffs === coeffs + # A refresh of the same shape is a write, not a rebuild: the interpolations read + # the panel's own storage, so the objects survive and nothing is allocated. + set_polar!(panel, alphas, [0.1, 0.4, 0.7, 1.0, 0.9], fill(0.02, 5), fill(-0.1, 5)) + cl_interp = panel.cl_interp + set_polar!(panel, alphas, [0.2, 0.5, 0.8, 1.1, 1.0], fill(0.03, 5), fill(-0.2, 5)) + @test panel.cl_interp === cl_interp + @test calculate_cl(panel, alphas[2]) ≈ 0.5 + cl_next, cd_next, cm_next = [0.2, 0.5, 0.8, 1.1, 1.0], fill(0.03, 5), fill(-0.2, 5) + refresh(p, a, l, d, m) = set_polar!(p, a, l, d, m) + refresh(panel, alphas, cl_next, cd_next, cm_next) + @test (@allocated refresh(panel, alphas, cl_next, cd_next, cm_next)) == 0 + @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)) From a8db009ef8c75cd65ab4317cb7039c1be04d5e0e Mon Sep 17 00:00:00 2001 From: Bart Date: Tue, 1 Sep 2026 10:16:28 +0200 Subject: [PATCH 2/5] Let a live polar refresh run without touching the heap A refresh over 60 panels allocated 15.3 MB, nearly all of it inside NeuralFoil's forward pass: every layer built a fresh matrix for its product and another for its activation, the Mahalanobis distance sliced a column and did a matvec per case, and the symmetry embedding ran the whole thing again for the flipped inputs before combining the two into a third matrix. NeuralFoilWorkspace holds one set of layer activations per symmetry plus the flipped input batch, and fused_output! runs both passes through it. Bias and activation fold into one in-place sweep, the penalty is subtracted case by case straight onto the confidence logit, and the flip is undone while averaging rather than into a third matrix. flipped_row is now the one place the output mirror symmetry is written down. LivePolars caches the model and that workspace and preallocates every per-panel buffer a refresh touches, so shapes deform in place, the coefficients decode in place and the panel takes them as views. The surface-pressure pass shares the same workspace and reconstructs Cp through reused buffers. set_polar! had two leaks that only showed under real use: its getfield-by-symbol loop went dynamic once the four sources were no longer the same type, and storing live_shape boxes the immutable on its way into a Union field, which cost 48 bytes a panel even when the shape had not changed. refresh_live_polars! now allocates nothing, refresh_live_pressure! 96% less, and polar_drift, live_surface_friction! and live_shape_offset! nothing. Run time is unchanged; the forward pass is what is left. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_0173NMo5G7SLMYjzCfr8L41u --- CHANGELOG.md | 24 ++ docs/make.jl | 2 +- docs/src/private_functions.md | 22 ++ docs/src/private_types.md | 2 + src/airfoil_aero/AirfoilAero.jl | 2 +- .../airfoil_solvers/neuralfoil_solver.jl | 171 +++++++-- src/airfoil_aero/deform.jl | 56 ++- src/airfoil_aero/live_polar.jl | 153 ++++++--- src/airfoil_aero/neuralfoil.jl | 325 ++++++++++++++---- src/panel.jl | 45 ++- test/airfoil_aero/test_live_polar.jl | 19 +- 11 files changed, 658 insertions(+), 163 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a072d783..ae33c142 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,30 @@ update: `set_polar!` on an unchanged table shape allocates nothing and leaves the interpolation objects in place. A table that changes shape is rebuilt once. +- `refresh_live_polars!` now allocates nothing at all, where a refresh over 60 + panels cost 15.3 MB before; the forward pass itself is what a refresh costs + now, and its run time is unchanged. NeuralFoil's pass was the bulk of the + garbage: every layer built a fresh matrix for its product, another for its + activation, and the symmetry embedding did the whole thing twice more for the + flipped case. `NeuralFoilWorkspace` holds one set of layer activations per + symmetry and `fused_output!` runs the pass through them, with the bias and + activation folded into one in-place sweep, the Mahalanobis penalty subtracted + from the confidence logit case by case rather than through a vector of + distances, and the flip undone while averaging instead of into a third matrix. + `LivePolars` keeps that workspace, the network itself, and a buffer for every + per-panel quantity a refresh touches, so the shapes deform in place + (`deform_kulfan!`), the coefficients decode in place (`decode_coefficients!`) + and the panel takes them as views. `polar_drift`, `live_surface_friction!` and + `live_shape_offset!` allocate nothing either. +- `refresh_live_pressure!` allocates 96% less (3.7 MB to 146 kB a refresh over + 60 panels): it shares the polar refresh's workspace, reads the edge velocities + out of it into storage it holds, and reconstructs each panel's `Cp` through + one set of reused buffers (`contour_arc!`, `contour_pressure!`). What is left + is the monotone interpolation object itself, one per panel. +- `set_polar!` no longer allocates when handed views rather than vectors, and + stores the `live_shape` only when the panel is not already carrying it — an + immutable that size is boxed on its way into a `Union` field, which cost a + refresh 48 bytes a panel even when the shape had not changed. ## VortexStepMethod v4.3.0 2026-08-31 diff --git a/docs/make.jl b/docs/make.jl index 2688d0c7..fd446529 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -14,7 +14,7 @@ makedocs(; VortexStepMethod.ObjAdapter], authors="Uwe Fechner , Bart van de Lint and contributors", sitename="VortexStepMethod.jl", - format = Documenter.HTML(prettyurls = haskey(ENV, "CI")), + format = Documenter.HTML(prettyurls = haskey(ENV, "CI"), size_threshold = 400_000), pages=[ "Home" => "index.md", "How it works" => "explanation.md", diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index fa41f393..197866a0 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -100,6 +100,7 @@ update_non_deformed_sections! calculate_new_aero_data set_polar! refresh_polar! +polar_column! polar_model reads_from same_knots @@ -134,6 +135,8 @@ CurrentModule = VortexStepMethod.AirfoilAero ### Kulfan CST parametrization ```@docs +deform_kulfan! +chord_residual! bernstein_basis class_function leading_edge_basis @@ -161,14 +164,27 @@ load_neuralfoil_model neuralfoil_section neuralfoil_fused_output fused_output +fused_output! decode_surface_velocity +decode_surface_velocity! +surface_velocity_rows decode_coefficients +decode_coefficients! nn_forward +nn_forward! +add_bias! +layer_buffers +case_capacity prepare_inputs fill_case_input! flip_inputs +flip_inputs! flip_outputs +flipped_row +fuse_flipped! squared_mahalanobis_distance +mahalanobis_case +penalize_confidence! swish sigmoid ``` @@ -187,10 +203,16 @@ write_node_table flat_plate_cf neuralfoil_contour_solution contour_arc +contour_arc! arc_at_chord +arc_at_chord! +sorted_surface! +sort_pairs! trailing_edge_speed velocity_knots +velocity_knots! contour_pressure +contour_pressure! fill_node_nans! live_xfoil_solver xfoil_ramp diff --git a/docs/src/private_types.md b/docs/src/private_types.md index ff99a2c5..0941c427 100644 --- a/docs/src/private_types.md +++ b/docs/src/private_types.md @@ -31,6 +31,8 @@ LivePolarSettings LivePolars NeuralFoilModel NeuralFoilResult +NeuralFoilWorkspace +ContourPressureScratch ``` ```@meta CurrentModule = VortexStepMethod diff --git a/src/airfoil_aero/AirfoilAero.jl b/src/airfoil_aero/AirfoilAero.jl index 140da02a..16f10813 100644 --- a/src/airfoil_aero/AirfoilAero.jl +++ b/src/airfoil_aero/AirfoilAero.jl @@ -16,11 +16,11 @@ 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") include("airfoil_solvers/neuralfoil_solver.jl") +include("live_polar.jl") include("airfoil_io.jl") include("polar_gen.jl") include("polar_export.jl") diff --git a/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl b/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl index f49eed1e..54f3ba1e 100644 --- a/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl +++ b/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl @@ -49,8 +49,17 @@ along the upper surface toward the trailing edge, negative along the lower. The natural coordinate across a blunt nose, where a small chordwise step is a long one along the skin. """ -function contour_arc(x, y, le) - arc = zeros(Float64, length(x)) +contour_arc(x, y, le) = contour_arc!(zeros(Float64, length(x)), x, y, le) + +""" + contour_arc!(arc, x, y, le) -> arc + +[`contour_arc`](@ref) written into `arc`, which is resized to the contour. A live +pressure refresh walks panel after panel through one buffer this way. +""" +function contour_arc!(arc::Vector{Float64}, x, y, le) + resize!(arc, length(x)) + arc[le] = 0.0 for k in (le - 1):-1:1 arc[k] = arc[k + 1] + hypot(x[k] - x[k + 1], y[k] - y[k + 1]) end @@ -60,6 +69,31 @@ function contour_arc(x, y, le) return arc end +""" + ContourPressureScratch() + +The buffers a surface-pressure reconstruction reuses, so [`contour_pressure!`](@ref) +allocates nothing per panel or per frame. Each grows to whatever contour and station +count it is first handed and stays that size. +""" +struct ContourPressureScratch + "One surface's node chord fractions, ascending." + surface_x::Vector{Float64} + "The same nodes' signed arc length." + surface_arc::Vector{Float64} + "Arc length of the network's stations on the upper surface." + upper_arc::Vector{Float64} + "Arc length of the network's stations on the lower surface." + lower_arc::Vector{Float64} + "Arc length of every knot of the joined edge-velocity curve." + knots::Vector{Float64} + "Signed edge velocity at those knots." + values::Vector{Float64} +end + +ContourPressureScratch() = ContourPressureScratch(Float64[], Float64[], Float64[], + Float64[], Float64[], Float64[]) + """ arc_at_chord(x, arc, indices, fractions) -> Vector @@ -67,14 +101,66 @@ Signed arc length at each chord fraction, read off the contour nodes `indices` one surface the fractions belong to. """ function arc_at_chord(x, arc, indices, fractions) - nodes = collect(indices) - order = sortperm(x[nodes]) - xs, as = x[nodes][order], arc[nodes][order] - return [begin + scratch = ContourPressureScratch() + return arc_at_chord!(zeros(length(fractions)), scratch, x, arc, indices, fractions) +end + +""" + arc_at_chord!(out, scratch, x, arc, indices, fractions) -> out + +[`arc_at_chord`](@ref) written into `out`, taking the sorted copy of the surface from +`scratch` rather than allocating one per call. +""" +function arc_at_chord!(out::Vector{Float64}, scratch::ContourPressureScratch, x, arc, + indices, fractions) + xs, as = sorted_surface!(scratch, x, arc, indices) + resize!(out, length(fractions)) + @inbounds for (i, f) in enumerate(fractions) j = clamp(searchsortedfirst(xs, f), 2, length(xs)) gap = max(xs[j] - xs[j - 1], eps()) - as[j - 1] + (f - xs[j - 1]) / gap * (as[j] - as[j - 1]) - end for f in fractions] + out[i] = as[j - 1] + (f - xs[j - 1]) / gap * (as[j] - as[j - 1]) + end + return out +end + +""" + sorted_surface!(scratch, x, arc, indices) -> (surface_x, surface_arc) + +One surface's nodes as chord fraction and signed arc length, ascending in chord, in the +scratch's own storage. Sorted by insertion, which a surface listed from one end to the +other is either already in or exactly reversed from, and which needs no scratch of its +own. +""" +function sorted_surface!(scratch::ContourPressureScratch, x, arc, indices) + surface_x, surface_arc = scratch.surface_x, scratch.surface_arc + n = length(indices) + resize!(surface_x, n) + resize!(surface_arc, n) + @inbounds for (k, node) in enumerate(indices) + surface_x[k] = x[node] + surface_arc[k] = arc[node] + end + sort_pairs!(surface_x, surface_arc) + return (surface_x, surface_arc) +end + +""" + sort_pairs!(sorted, carried) + +Insertion-sort `sorted` ascending, carrying `carried` along with it. Stable, in place, +and linear on input that is already ordered, which both of its callers hand it. +""" +function sort_pairs!(sorted::Vector{Float64}, carried::Vector{Float64}) + @inbounds for k in 2:length(sorted) + key, along = sorted[k], carried[k] + j = k - 1 + while j >= 1 && sorted[j] > key + sorted[j + 1], carried[j + 1] = sorted[j], carried[j] + j -= 1 + end + sorted[j + 1], carried[j + 1] = key, along + end + return nothing end """ @@ -111,20 +197,51 @@ than extrapolated off the end of one. The stagnation point enters as its own kno where `ue` interpolated linearly between those two stations reaches zero — the stagnation-point-flow result, `ue` being linear in arc length there. """ -function velocity_knots(station_x, ue_upper, ue_lower, x, arc, le) - upper_arc = arc_at_chord(x, arc, 1:le, station_x) - lower_arc = arc_at_chord(x, arc, le:length(x), station_x) +velocity_knots(station_x, ue_upper, ue_lower, x, arc, le) = + velocity_knots!(ContourPressureScratch(), station_x, ue_upper, ue_lower, x, arc, le) + +""" + velocity_knots!(scratch, station_x, ue_upper, ue_lower, x, arc, le) + -> (knots, values) + +[`velocity_knots`](@ref) assembled in `scratch`. The two vectors returned are the +scratch's own and are overwritten by the next call. +""" +function velocity_knots!(scratch::ContourPressureScratch, station_x, ue_upper, ue_lower, + x, arc, le) + n_stations = length(station_x) + upper_arc = arc_at_chord!(scratch.upper_arc, scratch, x, arc, 1:le, station_x) + lower_arc = arc_at_chord!(scratch.lower_arc, scratch, x, arc, le:length(x), + station_x) speed = trailing_edge_speed(ue_upper, ue_lower, upper_arc, lower_arc, arc[1], arc[end]) lower_first, upper_first = abs(ue_lower[1]), abs(ue_upper[1]) stagnation = lower_arc[1] + lower_first / max(lower_first + upper_first, eps()) * (upper_arc[1] - lower_arc[1]) - knots = [arc[end]; reverse(lower_arc); stagnation; upper_arc; arc[1]] - values = [-speed; reverse(-abs.(ue_lower)); 0.0; abs.(ue_upper); speed] - order = sortperm(knots) - knots, values = knots[order], values[order] - keep = [1; [k for k in 2:length(knots) if knots[k] > knots[k - 1] + 1e-9]] - return knots[keep], values[keep] + knots = resize!(scratch.knots, 2 * n_stations + 3) + values = resize!(scratch.values, 2 * n_stations + 3) + knots[1], values[1] = arc[end], -speed + for k in 1:n_stations + knots[1 + k], values[1 + k] = lower_arc[n_stations + 1 - k], + -abs(ue_lower[n_stations + 1 - k]) + end + knots[n_stations + 2], values[n_stations + 2] = stagnation, 0.0 + for k in 1:n_stations + knots[n_stations + 2 + k], values[n_stations + 2 + k] = upper_arc[k], + abs(ue_upper[k]) + end + knots[end], values[end] = arc[1], speed + sort_pairs!(knots, values) + kept, previous = 1, knots[1] + for k in 2:length(knots) + current = knots[k] + if current > previous + 1e-9 + kept += 1 + knots[kept], values[kept] = current, values[k] + end + previous = current + end + return (resize!(knots, kept), resize!(values, kept)) end """ @@ -135,10 +252,24 @@ end one-curve edge velocity ([`velocity_knots`](@ref)), interpolates it with a shape-preserving monotone cubic in arc length, and squares it into `Cp = 1 - ue²`. """ -function contour_pressure(station_x, ue_upper, ue_lower, x, arc, le) - knots, values = velocity_knots(station_x, ue_upper, ue_lower, x, arc, le) +contour_pressure(station_x, ue_upper, ue_lower, x, arc, le) = + contour_pressure!(zeros(length(x)), ContourPressureScratch(), station_x, ue_upper, + ue_lower, x, arc, le) + +""" + contour_pressure!(cp, scratch, station_x, ue_upper, ue_lower, x, arc, le) -> cp + +[`contour_pressure`](@ref) written into `cp`, with `scratch` carrying the edge-velocity +curve. Only the monotone interpolation itself still allocates, once per contour. +""" +function contour_pressure!(cp, scratch::ContourPressureScratch, station_x, ue_upper, + ue_lower, x, arc, le) + knots, values = velocity_knots!(scratch, station_x, ue_upper, ue_lower, x, arc, le) speed = interpolate(knots, values, FritschButlandMonotonicInterpolation()) - return [1 - speed(clamp(a, first(knots), last(knots)))^2 for a in arc] + @inbounds for k in eachindex(cp, arc) + cp[k] = 1 - speed(clamp(arc[k], first(knots), last(knots)))^2 + end + return cp end """ diff --git a/src/airfoil_aero/deform.jl b/src/airfoil_aero/deform.jl index a911a565..8bd97fe6 100644 --- a/src/airfoil_aero/deform.jl +++ b/src/airfoil_aero/deform.jl @@ -69,16 +69,46 @@ 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) + out = KulfanParameters(copy(base.upper_weights), copy(base.lower_weights), + base.leading_edge_weight, base.TE_thickness) + return deform_kulfan!(out, basis, base, upper_deflection, lower_deflection, + similar(basis.x)) +end + +""" + deform_kulfan!(out, basis, base, upper_deflection, lower_deflection, residual) + -> out + +[`deform_kulfan`](@ref) written into the weight vectors `out` already holds, with +`residual` as the scratch [`chord_residual!`](@ref) fills — the form a live polar +refreshes a panel's shape with, since it neither allocates nor changes the object the +panel points at. + +Only the weights are written. The leading-edge weight and trailing-edge thickness are +carried by a deformation unchanged, so `out` must already agree with `base` on both; +it does when `out` was built from `base`, which is checked. +""" +function deform_kulfan!(out::KulfanParameters, basis::KulfanBasis, + base::KulfanParameters, upper_deflection::AbstractVector, + lower_deflection::AbstractVector, residual::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 * chord_residual(basis, upper_deflection), - base.lower_weights .+ basis.projection * chord_residual(basis, lower_deflection), - base.leading_edge_weight, base.TE_thickness) + (out.leading_edge_weight == base.leading_edge_weight && + out.TE_thickness == base.TE_thickness) || throw(ArgumentError( + "deform_kulfan! writes weights only, so its target must carry the base " * + "airfoil's leading-edge weight and trailing-edge thickness.")) + for (weights, base_weights, deflection) in + ((out.upper_weights, base.upper_weights, upper_deflection), + (out.lower_weights, base.lower_weights, lower_deflection)) + chord_residual!(residual, basis, deflection) + mul!(weights, basis.projection, residual) + weights .+= base_weights + end + return out end """ @@ -94,9 +124,19 @@ magnitude past the ones it is correcting, and the airfoil that comes back is not Removes the straight line through the deflection's own endpoints, which [`chord_line`](@ref) returns, and gives back what is left. """ -function chord_residual(basis::KulfanBasis, deflection::AbstractVector) +chord_residual(basis::KulfanBasis, deflection::AbstractVector) = + chord_residual!(similar(basis.x), basis, deflection) + +""" + chord_residual!(residual, basis, deflection) -> residual + +[`chord_residual`](@ref) written into storage the caller owns. +""" +function chord_residual!(residual::AbstractVector, basis::KulfanBasis, + deflection::AbstractVector) offset, slope = chord_line(basis, deflection) - return deflection .- (offset .+ slope .* basis.x) + residual .= deflection .- (offset .+ slope .* basis.x) + return residual end """ @@ -119,6 +159,10 @@ end deform_kulfan(basis::KulfanBasis, base::KulfanParameters, camber::AbstractVector) = deform_kulfan(basis, base, camber, camber) +deform_kulfan!(out::KulfanParameters, basis::KulfanBasis, base::KulfanParameters, + camber::AbstractVector, residual::AbstractVector) = + deform_kulfan!(out, basis, base, camber, camber, residual) + """ control_point_deflection(basis, fractions, deflections) -> Vector{Float64} diff --git a/src/airfoil_aero/live_polar.jl b/src/airfoil_aero/live_polar.jl index 215eccbc..2d4e429f 100644 --- a/src/airfoil_aero/live_polar.jl +++ b/src/airfoil_aero/live_polar.jl @@ -32,9 +32,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 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 +Holds the fixed CST basis, the per-panel reference angles, the network and every buffer a +refresh needs, so a refresh writes through storage that is already there and allocates +nothing at all. 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 @@ -44,7 +45,7 @@ mutable struct LivePolars basis::KulfanBasis "Undeformed Kulfan parameters per panel." base::Vector{KulfanParameters} - "Deformed Kulfan parameters per panel, rewritten every refresh." + "Deformed Kulfan parameters per panel, rewritten in place every refresh." deformed::Vector{KulfanParameters} "Reference angle [rad] the last refresh sampled about, per panel." alpha_ref::Vector{Float64} @@ -56,6 +57,34 @@ mutable struct LivePolars pressure_inputs::Matrix{Float32} "Lowest confidence NeuralFoil reported over each panel's samples, last refresh." confidence::Vector{Float64} + "The network the settings name, held so a refresh never goes to the model cache." + model::NeuralFoilModel + "Forward-pass scratch, sized for the sampled batch and so wide enough for either pass." + work::NeuralFoilWorkspace + "Per-panel Reynolds number of the current pass." + reynolds::Vector{Float64} + "Per-panel angle [rad] the surface-pressure pass was evaluated at." + alpha_at::Vector{Float64} + "Lift coefficient of every sample of the last refresh." + sample_cl::Vector{Float64} + "Drag coefficient of every sample of the last refresh." + sample_cd::Vector{Float64} + "Moment coefficient of every sample of the last refresh." + sample_cm::Vector{Float64} + "Analysis confidence of every sample of the last refresh." + sample_confidence::Vector{Float64} + "Upper-surface edge velocity per station and panel, last surface-pressure pass." + ue_upper::Matrix{Float64} + "Lower-surface edge velocity per station and panel, last surface-pressure pass." + ue_lower::Matrix{Float64} + "Signed arc length of the contour one panel's pressure is being spread over." + arc::Vector{Float64} + "Buffers the pressure reconstruction walks from panel to panel through." + pressure_scratch::ContourPressureScratch + "The representable part of one panel's deflection." + residual::Vector{Float64} + "One panel's weight change off its base shape." + weight_delta::Vector{Float64} end function LivePolars(base::AbstractVector{KulfanParameters}; @@ -69,13 +98,23 @@ function LivePolars(base::AbstractVector{KulfanParameters}; 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), - zeros(length(offsets)), - zeros(Float32, 25, n_panels * length(offsets)), - zeros(Float32, 25, n_panels), zeros(n_panels)) + n_panels, n_samples = length(base), length(offsets) + n_weights = length(base[1].upper_weights) + n_cases = n_panels * n_samples + model = load_neuralfoil_model(settings.model_size; weights_dir=settings.weights_dir) + work = NeuralFoilWorkspace(model, n_cases) + n_network_stations = length(work.stations) + deformed = [KulfanParameters(copy(p.upper_weights), copy(p.lower_weights), + p.leading_edge_weight, p.TE_thickness) for p in base] + return LivePolars(settings, KulfanBasis(; n_stations, n_weights), + collect(base), deformed, zeros(n_panels), zeros(n_samples), + zeros(Float32, 25, n_cases), zeros(Float32, 25, n_panels), + zeros(n_panels), model, work, zeros(n_panels), zeros(n_panels), + zeros(n_cases), zeros(n_cases), zeros(n_cases), zeros(n_cases), + zeros(n_network_stations, n_panels), + zeros(n_network_stations, n_panels), + Float64[], ContourPressureScratch(), zeros(n_stations), + zeros(n_weights)) end """ @@ -114,21 +153,30 @@ function polar_drift(live::LivePolars, alpha::AbstractVector) "polar_drift: $(length(alpha)) angles for $(length(live.alpha_ref)) panels.")) offsets = live.settings.offsets reach = min(-first(offsets), last(offsets)) - return maximum(abs.(alpha .- live.alpha_ref)) / reach + return maximum(i -> abs(alpha[i] - live.alpha_ref[i]), + eachindex(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. +`live.deformed`, which is returned. `nothing` writes the base shapes back. The +deformation is an analytic perturbation of one fixed weight vector, so it never refits +and never inherits the non-uniqueness of a fit. + +Each entry keeps the object it already was, its weights overwritten, so a panel holding +one from an earlier refresh follows the current shape and no frame allocates a new one. """ 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]) + shape, base = live.deformed[i], live.base[i] + if isnothing(deflection) + shape.upper_weights .= base.upper_weights + shape.lower_weights .= base.lower_weights + else + deform_kulfan!(shape, live.basis, base, deflection[i], live.residual) + end end return live.deformed end @@ -147,7 +195,7 @@ function apply_live_shapes!(live::LivePolars, panels; deflection=nothing) "base airfoils.")) deform_live_shapes!(live, deflection) for (panel, shape) in zip(panels, live.deformed) - panel.live_shape = shape + panel.live_shape === shape || (panel.live_shape = shape) end return nothing end @@ -162,9 +210,10 @@ deform the base airfoil by `deflection` (a chord-normalized deflection on `live. 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 the three interpolation rebuilds `Interpolations` needs to take the -new values, which it copies rather than references. +The write is in place — same knot count, same vectors — and so is everything around it: +the shapes, the network inputs, both symmetries of the forward pass and the decoded +coefficients all land in storage `live` already holds, 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. @@ -182,9 +231,8 @@ function refresh_live_polars!(live::LivePolars, panels, alpha_ref, reynolds; 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 + live.alpha_ref .= alpha_ref + live.reynolds .= reynolds offsets = live.settings.offsets n_samples = length(offsets) @@ -192,23 +240,24 @@ function refresh_live_polars!(live::LivePolars, panels, alpha_ref, reynolds; for i in 1:n_panels for k in 1:n_samples fill_case_input!(live.inputs, (i - 1) * n_samples + k, live.deformed[i], - rad2deg(alpha_vec[i] + offsets[k]), re_vec[i], + rad2deg(live.alpha_ref[i] + offsets[k]), live.reynolds[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)) + decode_coefficients!(live.sample_cl, live.sample_cd, live.sample_cm, + live.sample_confidence, + fused_output!(live.work, live.inputs, live.model)) for i in 1:n_panels samples = ((i - 1) * n_samples + 1):(i * n_samples) - live.knots .= alpha_vec[i] .+ offsets - live.confidence[i] = minimum(view(confidence, samples)) - set_polar!(panels[i], live.knots, view(cl, samples), view(cd, samples), - view(cm, samples); shape=live.deformed[i]) + live.knots .= live.alpha_ref[i] .+ offsets + live.confidence[i] = minimum(view(live.sample_confidence, samples)) + set_polar!(panels[i], live.knots, view(live.sample_cl, samples), + view(live.sample_cd, samples), view(live.sample_cm, samples); + shape=live.deformed[i]) end - return minimum(confidence) + return minimum(live.sample_confidence) end """ @@ -293,7 +342,8 @@ panel forces follow the deformed shape; without this the pattern that spreads th 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 +One batched forward pass over all panels — through the same workspace the polar refresh +uses, whose inputs it has already consumed — 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 the contours the traction pattern is indexed on — @@ -308,22 +358,20 @@ function refresh_live_pressure!(cp, live::LivePolars, contour_x, contour_y, 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) + live.alpha_at .= alpha + live.reynolds .= 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) + rad2deg(live.alpha_at[i]), live.reynolds[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, ue_upper, ue_lower = decode_surface_velocity( - fused_output(live.pressure_inputs, model)) + decode_surface_velocity!(live.ue_upper, live.ue_lower, + fused_output!(live.work, live.pressure_inputs, live.model)) for i in 1:n_panels - arc = contour_arc(contour_x[i], contour_y[i], leading_edge[i]) - cp[i] .= contour_pressure(station_x, view(ue_upper, :, i), - view(ue_lower, :, i), contour_x[i], arc, - leading_edge[i]) + arc = contour_arc!(live.arc, contour_x[i], contour_y[i], leading_edge[i]) + contour_pressure!(cp[i], live.pressure_scratch, live.work.stations, + view(live.ue_upper, :, i), view(live.ue_lower, :, i), + contour_x[i], arc, leading_edge[i]) end return cp end @@ -342,11 +390,12 @@ 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]) + re = reynolds isa Number ? float(reynolds) : float(reynolds[i]) + nodes = contour_x[i] + @inbounds for k in eachindex(cf[i], nodes) + cf[i][k] = flat_plate_cf(clamp(nodes[k], 0.0, 1.0), re) + end end return cf end @@ -380,8 +429,8 @@ function live_shape_offset!(offset, live::LivePolars, shape) "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) + live.weight_delta .= live.deformed[i].upper_weights .- live.base[i].upper_weights + mul!(offset[i], shape[i], live.weight_delta) end return offset end diff --git a/src/airfoil_aero/neuralfoil.jl b/src/airfoil_aero/neuralfoil.jl index 080a3919..4a229bdb 100644 --- a/src/airfoil_aero/neuralfoil.jl +++ b/src/airfoil_aero/neuralfoil.jl @@ -121,6 +121,50 @@ function load_neuralfoil_model(model_size::String="xlarge"; weights_dir=nothing) return model end +""" + layer_buffers(model, n_cases) -> Vector{Matrix{Float32}} + +One activation matrix per network layer, `n_cases` wide: the storage a forward pass +writes through. [`NeuralFoilWorkspace`](@ref) holds two sets of them, one per symmetry. +""" +layer_buffers(model::NeuralFoilModel, n_cases::Int) = + [zeros(Float32, size(w, 1), n_cases) for w in model.weights] + +""" + NeuralFoilWorkspace(model, n_cases) + +Scratch a symmetry-fused forward pass runs inside, so a caller evaluating a batch of the +same width every solve allocates nothing at all. Sized once against `model` for at most +`n_cases` cases; [`fused_output!`](@ref) takes any batch up to that width, which is what +lets one workspace serve both a sampled polar batch and the narrower pressure batch. +""" +struct NeuralFoilWorkspace + "Layer activations of the direct pass." + direct::Vector{Matrix{Float32}} + "Layer activations of the top/bottom-flipped pass." + flipped::Vector{Matrix{Float32}} + "The flipped input batch." + inputs_flip::Matrix{Float32} + "One input case less the training mean." + centered::Vector{Float64} + "NeuralFoil's boundary-layer station x/c, the axis a surface reconstruction reads." + stations::Vector{Float64} +end + +function NeuralFoilWorkspace(model::NeuralFoilModel, n_cases::Int) + n_cases >= 1 || throw(ArgumentError( + "A NeuralFoilWorkspace needs at least one case; got $n_cases.")) + n_outputs = size(last(model.weights), 1) + return NeuralFoilWorkspace(layer_buffers(model, n_cases), + layer_buffers(model, n_cases), + zeros(Float32, model.n_inputs, n_cases), + zeros(model.n_inputs), + compute_optimal_x_points((n_outputs - 6) ÷ 6)) +end + +"How many cases a workspace was sized for." +case_capacity(work::NeuralFoilWorkspace) = size(work.inputs_flip, 2) + """ nn_forward(x::AbstractMatrix, model::NeuralFoilModel) @@ -133,36 +177,92 @@ Forward pass through the neural network. # Returns - Output matrix of shape (n_outputs, n_cases) """ -function nn_forward(x::AbstractMatrix{T}, model::NeuralFoilModel) where T - n_layers = length(model.weights) +nn_forward(x::AbstractMatrix, model::NeuralFoilModel) = + Matrix(nn_forward!(layer_buffers(model, size(x, 2)), x, model)) - for i in 1:n_layers - x = model.weights[i] * x .+ model.biases[i] - if i < n_layers # No activation on last layer - x = swish(x) - end +""" + nn_forward!(layers, x, model, n_cases=size(x, 2)) -> AbstractMatrix + +Forward pass writing each layer's activations into `layers`, returning a view of the last +one over the `n_cases` columns used. The matrix products go through BLAS in place, so a +pass over a batch the buffers were sized for allocates nothing. +""" +function nn_forward!(layers::Vector{Matrix{Float32}}, x::AbstractMatrix, + model::NeuralFoilModel, n_cases::Int=size(x, 2)) + n_layers = length(model.weights) + out = view(layers[1], :, 1:n_cases) + mul!(out, model.weights[1], view(x, :, 1:n_cases)) + add_bias!(out, model.biases[1], n_layers > 1) + for i in 2:n_layers + out = view(layers[i], :, 1:n_cases) + mul!(out, model.weights[i], view(layers[i - 1], :, 1:n_cases)) + add_bias!(out, model.biases[i], i < n_layers) end + return view(layers[n_layers], :, 1:n_cases) +end - return x +""" + add_bias!(out, bias, activate) -> out + +Add a layer's bias to its activations in place, passing them through [`swish`](@ref) +unless this is the output layer, which carries no activation. +""" +function add_bias!(out::AbstractMatrix, bias::AbstractVector, activate::Bool) + @inbounds for case in axes(out, 2), k in eachindex(bias) + value = out[k, case] + bias[k] + out[k, case] = activate ? swish(value) : value + end + return out end """ squared_mahalanobis_distance(x::AbstractMatrix, model::NeuralFoilModel) -Compute squared Mahalanobis distance from training distribution. +Compute squared Mahalanobis distance from training distribution, one per case. This is used to penalize predictions far from the training data. """ function squared_mahalanobis_distance(x::AbstractMatrix, model::NeuralFoilModel) - n_cases = size(x, 2) - result = zeros(n_cases) + centered = zeros(model.n_inputs) + return [mahalanobis_case(x, case, model, centered) for case in axes(x, 2)] +end - for i in 1:n_cases - x_minus_mean = x[:, i] .- model.mean_inputs - result[i] = dot(x_minus_mean, model.inv_cov_inputs * x_minus_mean) +""" + mahalanobis_case(x, case, model, centered) -> Float64 + +One case's squared Mahalanobis distance, taking `centered` as the scratch the centred +input column is written into rather than allocating one per case. +""" +function mahalanobis_case(x::AbstractMatrix, case::Int, model::NeuralFoilModel, + centered::Vector{Float64}) + inv_cov = model.inv_cov_inputs + @inbounds for i in eachindex(centered) + centered[i] = x[i, case] - model.mean_inputs[i] end + distance = 0.0 + @inbounds for j in eachindex(centered) + weighted = 0.0 + for i in eachindex(centered) + weighted += inv_cov[i, j] * centered[i] + end + distance += weighted * centered[j] + end + return distance +end + +""" + penalize_confidence!(y, x, model, centered) -> y - return result +Subtract every case's Mahalanobis penalty from the confidence logit, row 1 of the network +output. This is what makes a shape far from the training distribution report a low +confidence, and it is applied to each symmetry before the two are fused. +""" +function penalize_confidence!(y::AbstractMatrix, x::AbstractMatrix, + model::NeuralFoilModel, centered::Vector{Float64}) + @inbounds for case in axes(y, 2) + y[1, case] -= mahalanobis_case(x, case, model, centered) / (2 * model.n_inputs) + end + return y end """ @@ -242,23 +342,27 @@ end Flip inputs for symmetry embedding (swap upper/lower, negate alpha). """ -function flip_inputs(x::AbstractMatrix{T}) where T - x_flip = copy(x) - - # Swap upper and lower weights with sign flip - x_flip[1:8, :] .= -x[9:16, :] # upper <- -lower - x_flip[9:16, :] .= -x[1:8, :] # lower <- -upper - - # Flip LE weight - x_flip[17, :] .= -x[17, :] - - # Flip sin(2*alpha) (index 19) - x_flip[19, :] .= -x[19, :] +flip_inputs(x::AbstractMatrix) = flip_inputs!(similar(x), x) - # Swap transition locations - x_flip[24, :] .= x[25, :] # xtr_upper <- xtr_lower - x_flip[25, :] .= x[24, :] # xtr_lower <- xtr_upper +""" + flip_inputs!(x_flip, x) -> x_flip +[`flip_inputs`](@ref) written into storage the caller owns: the upper and lower weights +swap with a sign change, the leading-edge weight and `sin(2·alpha)` negate, and the two +transition locations swap. Everything else carries over. +""" +function flip_inputs!(x_flip::AbstractMatrix, x::AbstractMatrix) + x_flip .= x + @inbounds for case in axes(x, 2) + for i in 1:8 + x_flip[i, case] = -x[8 + i, case] + x_flip[8 + i, case] = -x[i, case] + end + x_flip[17, case] = -x[17, case] + x_flip[19, case] = -x[19, case] + x_flip[24, case] = x[25, case] + x_flip[25, case] = x[24, case] + end return x_flip end @@ -267,30 +371,58 @@ end Flip outputs back after evaluating with flipped inputs. """ -function flip_outputs(y::AbstractMatrix{T}) where T - y_flip = copy(y) - N = (size(y, 1) - 6) ÷ 6 - - # CL (row 2) and CM (row 4) flip sign - y_flip[2, :] .= -y[2, :] - y_flip[4, :] .= -y[4, :] - - # Transition locations swap (rows 5, 6) - y_flip[5, :] .= y[6, :] - y_flip[6, :] .= y[5, :] - - if N > 0 - # Swap upper/lower theta+H blocks (rows 7:6+2N ↔ 7+3N:6+5N) - y_flip[7:(6+2N), :] .= y[(7+3N):(6+5N), :] - y_flip[(7+3N):(6+5N), :] .= y[7:(6+2N), :] - # Swap upper/lower ue/vinf blocks with sign flip (velocity mirrors) - y_flip[(7+2N):(6+3N), :] .= -1 .* y[(7+5N):(6+6N), :] - y_flip[(7+5N):(6+6N), :] .= -1 .* y[(7+2N):(6+3N), :] +function flip_outputs(y::AbstractMatrix) + n_stations = (size(y, 1) - 6) ÷ 6 + y_flip = similar(y) + for row in axes(y, 1) + source, sign = flipped_row(row, n_stations) + @views y_flip[row, :] .= sign .* y[source, :] end - return y_flip end +""" + flipped_row(row, n_stations) -> (source, sign) + +Which row of a network output a top/bottom-flipped case reads `row` from, and with which +sign. The single place the mirror symmetry of the output layout is written down, shared +by [`flip_outputs`](@ref) and [`fuse_flipped!`](@ref). + +`CL` and `CM` change sign, the two transition locations swap, and of the six +`n_stations`-long boundary-layer blocks the upper and lower halves swap — the momentum +thickness and shape factor as they are, the edge velocity mirrored, so it changes sign. +""" +function flipped_row(row::Int, n_stations::Int) + (row == 2 || row == 4) && return (row, -1) + row == 5 && return (6, 1) + row == 6 && return (5, 1) + (row <= 6 || n_stations == 0) && return (row, 1) + block = (row - 7) ÷ n_stations + block in (0, 1) && return (row + 3n_stations, 1) + block == 2 && return (row + 3n_stations, -1) + block in (3, 4) && return (row - 3n_stations, 1) + block == 5 && return (row - 3n_stations, -1) + return (row, 1) +end + +""" + fuse_flipped!(y, y_flip) -> y + +Average a direct output with its flipped counterpart in place, undoing the flip on the +way ([`flipped_row`](@ref)). Reading `y_flip` while writing `y` is what lets the fusion +land in storage that is already there instead of a third matrix. +""" +function fuse_flipped!(y::AbstractMatrix, y_flip::AbstractMatrix) + n_stations = (size(y, 1) - 6) ÷ 6 + @inbounds for row in axes(y, 1) + source, sign = flipped_row(row, n_stations) + for case in axes(y, 2) + y[row, case] = (y[row, case] + sign * y_flip[source, case]) / 2 + end + end + return y +end + """ compute_optimal_x_points(n) -> Vector{Float64} @@ -352,10 +484,25 @@ Turn a fused network output matrix into the integrated coefficients per case, un 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, :]))) + n_cases = size(y, 2) + return decode_coefficients!(zeros(n_cases), zeros(n_cases), zeros(n_cases), + zeros(n_cases), y) +end + +""" + decode_coefficients!(cl, cd, cm, confidence, y) -> (cl, cd, cm, confidence) + +[`decode_coefficients`](@ref) into vectors the caller owns, so a live polar reads a +refresh out of the network without allocating four vectors for it. +""" +function decode_coefficients!(cl, cd, cm, confidence, y::AbstractMatrix) + @inbounds for case in axes(y, 2) + cl[case] = y[2, case] / 2 + cd[case] = clamp(exp((y[3, case] - 2) * 2), 0.0, 1.0) + cm[case] = y[4, case] / 20 + confidence[case] = sigmoid(y[1, case]) + end + return (cl, cd, cm, confidence) end """ @@ -381,15 +528,29 @@ 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) +fused_output(x::AbstractMatrix, model::NeuralFoilModel) = + Matrix(fused_output!(NeuralFoilWorkspace(model, size(x, 2)), x, model)) - x_flip = flip_inputs(x) - y_flip = nn_forward(x_flip, model) - y_flip[1, :] .-= squared_mahalanobis_distance(x_flip, model) ./ (2 * model.n_inputs) +""" + fused_output!(work, x, model) -> AbstractMatrix - return (y .+ flip_outputs(y_flip)) ./ 2 +[`fused_output`](@ref) run entirely inside `work`, allocating nothing for a batch the +workspace has room for. The result is a view of the workspace's own storage and stays +valid until the next pass through it. +""" +function fused_output!(work::NeuralFoilWorkspace, x::AbstractMatrix, + model::NeuralFoilModel) + n_cases = size(x, 2) + n_cases <= case_capacity(work) || throw(ArgumentError( + "A NeuralFoilWorkspace sized for $(case_capacity(work)) cases was handed " * + "$n_cases.")) + x_flip = view(work.inputs_flip, :, 1:n_cases) + flip_inputs!(x_flip, x) + y = nn_forward!(work.direct, x, model, n_cases) + penalize_confidence!(y, x, model, work.centered) + y_flip = nn_forward!(work.flipped, x_flip, model, n_cases) + penalize_confidence!(y_flip, x_flip, model, work.centered) + return fuse_flipped!(y, y_flip) end """ @@ -431,10 +592,44 @@ fixed stations. Each `ue` matrix is `N × n_cases`. Interpolate these rather tha `Cp` is quadratic. """ function decode_surface_velocity(y::AbstractMatrix) - N = (size(y, 1) - 6) ÷ 6 - return (compute_optimal_x_points(N), - Matrix{Float64}(y[(7 + 2N):(6 + 3N), :]), - Matrix{Float64}(y[(7 + 5N):(6 + 6N), :])) + n_stations, upper, lower = surface_velocity_rows(y) + return (compute_optimal_x_points(n_stations), + Matrix{Float64}(y[upper, :]), Matrix{Float64}(y[lower, :])) +end + +""" + decode_surface_velocity!(ue_upper, ue_lower, y) -> (ue_upper, ue_lower) + +[`decode_surface_velocity`](@ref) into matrices the caller owns. The station axis is +fixed by the network and is carried by [`NeuralFoilWorkspace`](@ref) instead, so a live +pressure refresh reads a pass out without allocating anything for it. +""" +function decode_surface_velocity!(ue_upper::AbstractMatrix, ue_lower::AbstractMatrix, + y::AbstractMatrix) + _, upper, lower = surface_velocity_rows(y) + @inbounds for case in axes(y, 2) + for (k, row) in enumerate(upper) + ue_upper[k, case] = y[row, case] + end + for (k, row) in enumerate(lower) + ue_lower[k, case] = y[row, case] + end + end + return (ue_upper, ue_lower) +end + +""" + surface_velocity_rows(y) -> (n_stations, upper, lower) + +The number of boundary-layer stations a fused output matrix carries and the row ranges +its upper and lower edge-velocity ratios sit in. Viewing those rows is how a caller that +already holds the output reads the velocities without copying them out +([`decode_surface_velocity`](@ref) is the copying form). +""" +function surface_velocity_rows(y::AbstractMatrix) + n_stations = (size(y, 1) - 6) ÷ 6 + return (n_stations, (7 + 2n_stations):(6 + 3n_stations), + (7 + 5n_stations):(6 + 6n_stations)) end """ diff --git a/src/panel.jl b/src/panel.jl index 2c4de856..d703c4ce 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -297,19 +297,22 @@ Rewrite a panel's polar table with values at `alphas` [rad], ascending, and rebu interpolations. The table keeps the shape the panel was built with: a `POLAR_VECTORS` panel takes the angles as they are, and a `POLAR_MATRICES` panel takes them at every `delta` it already spans, since a regenerated polar carries its deflection as shape -rather than as a flap angle and so says the same thing at each. 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)). +rather than as a flap angle and so says the same thing at each. 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 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. +unchanged ([`polar_column!`](@ref)), and the interpolations read that storage rather than +a copy of it, so writing it is the whole update and the interpolation objects stay put +([`refresh_polar!`](@ref)). Only a table that changes shape is rebuilt. Values may be +views into a caller's own buffers; nothing here holds on to them. `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. +cannot drift apart. It is only stored when it is not already there: an immutable this +size is boxed on its way into a `Union` field, so the guard is what keeps a refresh of a +shape whose weights were rewritten in place free of that box. """ function set_polar!(panel::Panel, alphas, cl, cd, cm; shape=nothing) length(alphas) == length(cl) == length(cd) == length(cm) || @@ -323,16 +326,11 @@ function set_polar!(panel::Panel, alphas, cl, cd, cm; shape=nothing) panel.aero_model = polar_model(panel.cl_interp) 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)) - dst = getfield(panel, dst_sym) - if length(dst) == length(src) - dst .= src - else - setfield!(panel, dst_sym, collect(Float64, src)) - end - end + panel.live_shape === shape || (panel.live_shape = shape) + panel.alpha_knots = polar_column!(panel.alpha_knots, alphas) + panel.cl_coeffs = polar_column!(panel.cl_coeffs, cl) + panel.cd_coeffs = polar_column!(panel.cd_coeffs, cd) + panel.cm_coeffs = polar_column!(panel.cm_coeffs, cm) cl_interp = refresh_polar!(panel.cl_interp, panel.alpha_knots, panel.cl_coeffs) cd_interp = refresh_polar!(panel.cd_interp, panel.alpha_knots, panel.cd_coeffs) cm_interp = refresh_polar!(panel.cm_interp, panel.alpha_knots, panel.cm_coeffs) @@ -342,6 +340,19 @@ function set_polar!(panel::Panel, alphas, cl, cd, cm; shape=nothing) return nothing end +""" + polar_column!(stored, values) -> Vector{Float64} + +One column of a rewritten polar table in the panel's own `stored` vector, which is +returned, or in a fresh vector when the sample count has changed — the one case a +rewritten table has to grow, and the one that costs the panel's interpolations a rebuild. +""" +function polar_column!(stored::Vector{Float64}, values) + length(stored) == length(values) || return collect(Float64, values) + stored .= values + return stored +end + """ refresh_polar!(old, knots, values) -> Extrapolation diff --git a/test/airfoil_aero/test_live_polar.jl b/test/airfoil_aero/test_live_polar.jl index 61c9e869..43a21f02 100644 --- a/test/airfoil_aero/test_live_polar.jl +++ b/test/airfoil_aero/test_live_polar.jl @@ -226,6 +226,19 @@ end @test panels[1].live_shape === live.deformed[1] @test panels[1].live_shape.upper_weights ≈ live.base[1].upper_weights + # Which makes a refresh the forward pass and nothing else. The shapes, the network + # inputs, both symmetries of the pass and the decoded coefficients all land in + # storage `LivePolars` already holds, so a solve loop refreshing every step adds + # nothing to the heap — the property the whole struct is shaped around. + deflections = fill(camber, n_panels) + refresh_deformed(l, p, d) = refresh_live_polars!(l, p, deg2rad(6.0), 3e6; + deflection=d) + refresh_base(l, p) = refresh_live_polars!(l, p, deg2rad(6.0), 3e6) + refresh_deformed(live, panels, deflections) + @test (@allocated refresh_deformed(live, panels, deflections)) == 0 + refresh_base(live, panels) + @test (@allocated refresh_base(live, panels)) == 0 + @test_throws ArgumentError refresh_live_polars!(live, panels[1:2], 0.0, 3e6) @test_throws ArgumentError LivePolars(fill(base, 2); settings=LivePolarSettings(; offsets=deg2rad.([1.0, 2.0, 3.0]))) @@ -323,7 +336,11 @@ end 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]) + contours, reynolds = [contour_x, contour_x], [3e6, 1e6] + live_surface_friction!(cf, contours, reynolds) + friction(c, x, re) = live_surface_friction!(c, x, re) + friction(cf, contours, reynolds) + @test (@allocated friction(cf, contours, reynolds)) == 0 @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. From 4d7f9856f8ca853cfae4075c94aaafb81449d607 Mon Sep 17 00:00:00 2001 From: Bart Date: Tue, 1 Sep 2026 10:43:06 +0200 Subject: [PATCH 3/5] Make KulfanParameters mutable so a panel holds the shape, not a copy As an immutable it was boxed on its way into the panel's live_shape Union field, so storing it cost 48 bytes a panel per refresh even when the shape had not changed. That needed an === guard to avoid, and the guard only worked because === on an immutable is structural: the panel was really holding a copy that compared equal to the shape it was sampled from. Mutable, the field holds the object itself. The guard goes, deform_kulfan! writes all four fields rather than requiring its target to already agree with the base on two of them, and a source that rewrites a shape in place has by that act updated every panel flying it. Two separately built KulfanParameters no longer compare === or == on identical contents. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_0173NMo5G7SLMYjzCfr8L41u --- CHANGELOG.md | 13 +++++++++---- src/airfoil_aero/deform.jl | 17 +++++------------ src/airfoil_aero/live_polar.jl | 4 +++- src/panel.jl | 7 +++---- src/section_aero.jl | 7 ++++++- 5 files changed, 26 insertions(+), 22 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ae33c142..61bf21e8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -31,10 +31,15 @@ out of it into storage it holds, and reconstructs each panel's `Cp` through one set of reused buffers (`contour_arc!`, `contour_pressure!`). What is left is the monotone interpolation object itself, one per panel. -- `set_polar!` no longer allocates when handed views rather than vectors, and - stores the `live_shape` only when the panel is not already carrying it — an - immutable that size is boxed on its way into a `Union` field, which cost a - refresh 48 bytes a panel even when the shape had not changed. +- `set_polar!` no longer allocates when handed views rather than vectors. +- `KulfanParameters` is now mutable, so a live polar source rewrites one shape + per panel every solve instead of building a new one, and the panel pointing at + it follows without being told. As an immutable it was boxed on its way into + the panel's `live_shape` field, which cost a refresh 48 bytes a panel even + when the shape had not changed; a panel's `live_shape` is now the very object + it was sampled from rather than a copy that compares equal to it. Two + separately built `KulfanParameters` therefore no longer compare `===` or `==` + on identical contents. ## VortexStepMethod v4.3.0 2026-08-31 diff --git a/src/airfoil_aero/deform.jl b/src/airfoil_aero/deform.jl index 8bd97fe6..30814056 100644 --- a/src/airfoil_aero/deform.jl +++ b/src/airfoil_aero/deform.jl @@ -79,14 +79,9 @@ end deform_kulfan!(out, basis, base, upper_deflection, lower_deflection, residual) -> out -[`deform_kulfan`](@ref) written into the weight vectors `out` already holds, with -`residual` as the scratch [`chord_residual!`](@ref) fills — the form a live polar -refreshes a panel's shape with, since it neither allocates nor changes the object the -panel points at. - -Only the weights are written. The leading-edge weight and trailing-edge thickness are -carried by a deformation unchanged, so `out` must already agree with `base` on both; -it does when `out` was built from `base`, which is checked. +[`deform_kulfan`](@ref) written into the shape `out` already is, with `residual` as the +scratch [`chord_residual!`](@ref) fills — the form a live polar refreshes a panel's shape +with, since it neither allocates nor replaces the object the panel points at. """ function deform_kulfan!(out::KulfanParameters, basis::KulfanBasis, base::KulfanParameters, upper_deflection::AbstractVector, @@ -97,10 +92,8 @@ function deform_kulfan!(out::KulfanParameters, basis::KulfanBasis, (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.")) - (out.leading_edge_weight == base.leading_edge_weight && - out.TE_thickness == base.TE_thickness) || throw(ArgumentError( - "deform_kulfan! writes weights only, so its target must carry the base " * - "airfoil's leading-edge weight and trailing-edge thickness.")) + out.leading_edge_weight = base.leading_edge_weight + out.TE_thickness = base.TE_thickness for (weights, base_weights, deflection) in ((out.upper_weights, base.upper_weights, upper_deflection), (out.lower_weights, base.lower_weights, lower_deflection)) diff --git a/src/airfoil_aero/live_polar.jl b/src/airfoil_aero/live_polar.jl index 2d4e429f..ee98551b 100644 --- a/src/airfoil_aero/live_polar.jl +++ b/src/airfoil_aero/live_polar.jl @@ -174,6 +174,8 @@ function deform_live_shapes!(live::LivePolars, deflection) if isnothing(deflection) shape.upper_weights .= base.upper_weights shape.lower_weights .= base.lower_weights + shape.leading_edge_weight = base.leading_edge_weight + shape.TE_thickness = base.TE_thickness else deform_kulfan!(shape, live.basis, base, deflection[i], live.residual) end @@ -195,7 +197,7 @@ function apply_live_shapes!(live::LivePolars, panels; deflection=nothing) "base airfoils.")) deform_live_shapes!(live, deflection) for (panel, shape) in zip(panels, live.deformed) - panel.live_shape === shape || (panel.live_shape = shape) + panel.live_shape = shape end return nothing end diff --git a/src/panel.jl b/src/panel.jl index d703c4ce..8ddcd826 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -310,9 +310,8 @@ views into a caller's own buffers; nothing here holds on to them. `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. It is only stored when it is not already there: an immutable this -size is boxed on its way into a `Union` field, so the guard is what keeps a refresh of a -shape whose weights were rewritten in place free of that box. +cannot drift apart. The panel holds the shape itself, not a copy, so a source that +rewrites one in place has already updated every panel flying it. """ function set_polar!(panel::Panel, alphas, cl, cd, cm; shape=nothing) length(alphas) == length(cl) == length(cd) == length(cm) || @@ -326,7 +325,7 @@ function set_polar!(panel::Panel, alphas, cl, cd, cm; shape=nothing) panel.aero_model = polar_model(panel.cl_interp) panel.alpha_ref = (alphas[1] + alphas[end]) / 2 panel.alpha_window = (alphas[end] - alphas[1]) / 2 - panel.live_shape === shape || (panel.live_shape = shape) + panel.live_shape = shape panel.alpha_knots = polar_column!(panel.alpha_knots, alphas) panel.cl_coeffs = polar_column!(panel.cl_coeffs, cl) panel.cd_coeffs = polar_column!(panel.cd_coeffs, cd) diff --git a/src/section_aero.jl b/src/section_aero.jl index 30dee3fd..eb84ebb0 100644 --- a/src/section_aero.jl +++ b/src/section_aero.jl @@ -9,13 +9,18 @@ plot or a traction pattern reads without knowing how the shape was produced. Fit with `AirfoilAero.fit_kulfan_parameters`, deform one with `AirfoilAero.deform_kulfan` and turn one into coordinates with `AirfoilAero.kulfan_to_coordinates`. +Mutable, so a live polar source can rewrite one shape every solve instead of building a +new one, and the panel pointing at it follows without being told +(`AirfoilAero.deform_kulfan!`). That also makes a panel's `live_shape` the very object it +was sampled from rather than a copy that compares equal to it. + # 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 +mutable struct KulfanParameters upper_weights::Vector{Float64} lower_weights::Vector{Float64} leading_edge_weight::Float64 From 3e2d834aeff59dacad93917a4649ff83ef695a6a Mon Sep 17 00:00:00 2001 From: Bart Date: Tue, 1 Sep 2026 10:49:16 +0200 Subject: [PATCH 4/5] Note what the allocation work buys a multi-core sweep Single-core the refresh is arithmetic-bound and unchanged, so the changelog entry read as a wash. It is not: allocation is shared state and Julia's collector stops every worker, so a sweep running one model per core was losing half its wall time to GC. Eight concurrent workers go from 4.49 to 2.02 ms a refresh. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_0173NMo5G7SLMYjzCfr8L41u --- CHANGELOG.md | 31 +++++++++++++++++-------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 61bf21e8..bd09addb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,20 +12,23 @@ the interpolation objects in place. A table that changes shape is rebuilt once. - `refresh_live_polars!` now allocates nothing at all, where a refresh over 60 - panels cost 15.3 MB before; the forward pass itself is what a refresh costs - now, and its run time is unchanged. NeuralFoil's pass was the bulk of the - garbage: every layer built a fresh matrix for its product, another for its - activation, and the symmetry embedding did the whole thing twice more for the - flipped case. `NeuralFoilWorkspace` holds one set of layer activations per - symmetry and `fused_output!` runs the pass through them, with the bias and - activation folded into one in-place sweep, the Mahalanobis penalty subtracted - from the confidence logit case by case rather than through a vector of - distances, and the flip undone while averaging instead of into a third matrix. - `LivePolars` keeps that workspace, the network itself, and a buffer for every - per-panel quantity a refresh touches, so the shapes deform in place - (`deform_kulfan!`), the coefficients decode in place (`decode_coefficients!`) - and the panel takes them as views. `polar_drift`, `live_surface_friction!` and - `live_shape_offset!` allocate nothing either. + panels cost 15.3 MB before. On one core that buys little — the pass is + arithmetic-bound, not allocation-bound — but a sweep running one model per + core was spending half its wall time in the garbage collector, which stops + every worker. Eight workers refreshing concurrently went from 4.49 to 2.02 ms + a refresh, and the collector from 51% of the run to nothing. NeuralFoil's + forward pass was the bulk of the garbage: every layer built a fresh matrix for + its product, another for its activation, and the symmetry embedding did the + whole thing twice more for the flipped case. `NeuralFoilWorkspace` holds one + set of layer activations per symmetry and `fused_output!` runs the pass + through them, with the bias and activation folded into one in-place sweep, the + Mahalanobis penalty subtracted from the confidence logit case by case rather + than through a vector of distances, and the flip undone while averaging + instead of into a third matrix. `LivePolars` keeps that workspace, the network + itself, and a buffer for every per-panel quantity a refresh touches, so the + shapes deform in place (`deform_kulfan!`), the coefficients decode in place + (`decode_coefficients!`) and the panel takes them as views. `polar_drift`, + `live_surface_friction!` and `live_shape_offset!` allocate nothing either. - `refresh_live_pressure!` allocates 96% less (3.7 MB to 146 kB a refresh over 60 panels): it shares the polar refresh's workspace, reads the edge velocities out of it into storage it holds, and reconstructs each panel's `Cp` through From 1f79c948260239d7a54de94b369228a05eaef884 Mon Sep 17 00:00:00 2001 From: Bart Date: Tue, 1 Sep 2026 11:34:21 +0200 Subject: [PATCH 5/5] Measure the refresh allocations without boxing its return value On Julia 1.11 @allocated boxes a Float64 result, so the two new assertions read 16 bytes that refresh_live_polars! never allocated and failed CI there; 1.12 elides the box, which is why they passed locally. Reproduced on 1.11.9: measured returning the confidence, 16 bytes at top level and inside a let alike; measured returning nothing, 0 both ways. The closure had nothing to do with it. The two helpers now drop the confidence before the measurement, which is also why the set_polar! and live_surface_friction! assertions were green -- they measure calls returning nothing and a Vector. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_0173NMo5G7SLMYjzCfr8L41u --- test/airfoil_aero/test_live_polar.jl | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/test/airfoil_aero/test_live_polar.jl b/test/airfoil_aero/test_live_polar.jl index 43a21f02..9a799ec4 100644 --- a/test/airfoil_aero/test_live_polar.jl +++ b/test/airfoil_aero/test_live_polar.jl @@ -19,6 +19,22 @@ function polar_panels(n_panels; n_samples=5) end end +"""One `refresh_live_polars!` at the reference condition, its confidence dropped.""" +function refresh_once!(live, panels, deflection) + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6; deflection) + return nothing +end + +""" +Allocations of one `refresh_once!`, warmed up first. The confidence the refresh returns +has to be dropped rather than measured: on Julia 1.11 `@allocated` boxes a `Float64` +result, which reads as 16 bytes the refresh itself never allocated. +""" +function refresh_allocs(live, panels; deflection=nothing) + refresh_once!(live, panels, deflection) + return @allocated refresh_once!(live, panels, deflection) +end + @testset "Kulfan deformation" begin basis = KulfanBasis() base = KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) @@ -230,14 +246,8 @@ end # inputs, both symmetries of the pass and the decoded coefficients all land in # storage `LivePolars` already holds, so a solve loop refreshing every step adds # nothing to the heap — the property the whole struct is shaped around. - deflections = fill(camber, n_panels) - refresh_deformed(l, p, d) = refresh_live_polars!(l, p, deg2rad(6.0), 3e6; - deflection=d) - refresh_base(l, p) = refresh_live_polars!(l, p, deg2rad(6.0), 3e6) - refresh_deformed(live, panels, deflections) - @test (@allocated refresh_deformed(live, panels, deflections)) == 0 - refresh_base(live, panels) - @test (@allocated refresh_base(live, panels)) == 0 + @test refresh_allocs(live, panels; deflection=fill(camber, n_panels)) == 0 + @test refresh_allocs(live, panels) == 0 @test_throws ArgumentError refresh_live_polars!(live, panels[1:2], 0.0, 3e6) @test_throws ArgumentError LivePolars(fill(base, 2);