From a81fbe7a9029b5e0c3b5d3dec4b889839da2052b Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Tue, 22 Sep 2026 00:28:09 +0200 Subject: [PATCH 1/6] Check POLAR_MATRICES forwarddiff on affine tables instead of NeuralFoil's The ram-air wing's tip sections are ~70%-thick slices (#361) on which NeuralFoil returns Cd=1 at every alpha, and those outputs differ with the runner's CPU, so the LOOP solve lands on one of two different solutions depending on the runner. Tables that are affine in (alpha, delta) make the bilinear interpolation exact, keep every panel attached, and are the same bits on every machine. Co-Authored-By: Claude Opus 5 --- test/solver/test_forwarddiff.jl | 40 ++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 16 deletions(-) diff --git a/test/solver/test_forwarddiff.jl b/test/solver/test_forwarddiff.jl index 62792c2f..25f1ba01 100644 --- a/test/solver/test_forwarddiff.jl +++ b/test/solver/test_forwarddiff.jl @@ -5,6 +5,25 @@ using Test relative_error(jac, reference) = maximum(abs.(jac .- reference)) / maximum(abs, reference) +# Four-section arc whose POLAR_MATRICES tables are affine in (alpha, delta), so the +# bilinear interpolation is exact and has no kinks at the grid knots. +function affine_matrix_wing(n_panels) + alphas = deg2rad.(-5:5:25) + deltas = deg2rad.(-3:3:3) + cl = [0.2 + 5.5alpha + 1.5delta for alpha in alphas, delta in deltas] + cd = [0.03 + 0.2alpha + 0.05delta for alpha in alphas, delta in deltas] + cm = [-0.05 - 0.1alpha - 0.3delta for alpha in alphas, delta in deltas] + wing = Wing(n_panels, spanwise_distribution=LINEAR) + radius = 3.0 + for phi in deg2rad.((50, 17, -17, -50)) + le = [0.0, radius * sin(phi), radius * (cos(phi) - 1)] + add_section!(wing, le, le .+ [1.0, 0.0, 0.0], POLAR_MATRICES, + (alphas, deltas, cl, cd, cm)) + end + refine!(wing) + return wing +end + @testset "ForwardDiff linearize" begin n_panels = 10 span = 20.0 @@ -57,20 +76,9 @@ relative_error(jac, reference) = maximum(abs.(jac .- reference)) / maximum(abs, end @testset "AutoForwardDiff matches AutoFiniteDiff (LOOP, POLAR_MATRICES)" begin - # At this operating point the LOOP solve converges the mid-span - # panels' local alpha to ~6.5-9.1deg and the tip panels' to - # ~14.9-15.0deg. The default 5deg-spaced alpha_range=-5:5:15 puts the - # tip cluster within ulp-scale distance of the 15deg knot on some - # platforms (Windows, Julia 1.12), so AutoForwardDiff and the FD - # reference land on opposite sides of that non-differentiable kink - # and disagree by several percent (#360). Offset the grid so no knot - # is near either cluster (margin >0.8deg here vs <0.04deg before). - ram_wing = ram_air_matrix_wing(; n_panels=8, n_sections=4, - alpha_range=deg2rad.(-2:5:23), - delta_range=deg2rad.(-3:3:3), - ) - ram_body = BodyAerodynamics([ram_wing]) - ram_solver = Solver(ram_wing.n_panels, ram_wing.n_unrefined_sections; + matrix_wing = affine_matrix_wing(8) + matrix_body = BodyAerodynamics([matrix_wing]) + matrix_solver = Solver(matrix_wing.n_panels, matrix_wing.n_unrefined_sections; aerodynamic_model_type=VSM, rtol=1e-11, solver_type=LOOP, @@ -84,13 +92,13 @@ relative_error(jac, reference) = maximum(abs.(jac .- reference)) / maximum(abs, zeros(3)] jac_fwd, _, conv_fwd = VortexStepMethod.linearize( - ram_solver, ram_body, y_op; + matrix_solver, matrix_body, y_op; theta_idxs=1:4, va_vec_idxs=5:7, omega_idxs=8:10, aero_coeffs=true, backend=AutoForwardDiff()) @test conv_fwd jac_fd, _, conv_fd = VortexStepMethod.linearize( - ram_solver, ram_body, y_op; + matrix_solver, matrix_body, y_op; theta_idxs=1:4, va_vec_idxs=5:7, omega_idxs=8:10, aero_coeffs=true, backend=nothing) @test conv_fd From 190694dd6ba135d475f9cb0d07a3486b4f7fb69b Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Tue, 22 Sep 2026 01:06:44 +0200 Subject: [PATCH 2/6] Keep affine_matrix_wing's comment to one line and its panel count inside it Co-Authored-By: Claude Opus 5 --- test/solver/test_forwarddiff.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/test/solver/test_forwarddiff.jl b/test/solver/test_forwarddiff.jl index 25f1ba01..6fc3d4f6 100644 --- a/test/solver/test_forwarddiff.jl +++ b/test/solver/test_forwarddiff.jl @@ -5,15 +5,14 @@ using Test relative_error(jac, reference) = maximum(abs.(jac .- reference)) / maximum(abs, reference) -# Four-section arc whose POLAR_MATRICES tables are affine in (alpha, delta), so the -# bilinear interpolation is exact and has no kinks at the grid knots. -function affine_matrix_wing(n_panels) +# Affine tables: bilinear interpolation is exact, so the grid knots are not kinks. +function affine_matrix_wing() alphas = deg2rad.(-5:5:25) deltas = deg2rad.(-3:3:3) cl = [0.2 + 5.5alpha + 1.5delta for alpha in alphas, delta in deltas] cd = [0.03 + 0.2alpha + 0.05delta for alpha in alphas, delta in deltas] cm = [-0.05 - 0.1alpha - 0.3delta for alpha in alphas, delta in deltas] - wing = Wing(n_panels, spanwise_distribution=LINEAR) + wing = Wing(8, spanwise_distribution=LINEAR) radius = 3.0 for phi in deg2rad.((50, 17, -17, -50)) le = [0.0, radius * sin(phi), radius * (cos(phi) - 1)] @@ -76,7 +75,7 @@ end end @testset "AutoForwardDiff matches AutoFiniteDiff (LOOP, POLAR_MATRICES)" begin - matrix_wing = affine_matrix_wing(8) + matrix_wing = affine_matrix_wing() matrix_body = BodyAerodynamics([matrix_wing]) matrix_solver = Solver(matrix_wing.n_panels, matrix_wing.n_unrefined_sections; aerodynamic_model_type=VSM, From 06f88f3df17c658db4a0656406cf2cb0d313e99f Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Tue, 22 Sep 2026 14:49:44 +0200 Subject: [PATCH 3/6] Keep calc_forces!' per-panel loads out of line, which Julia 1.12 miscompiles for Dual{10} on Zen 4/5 panel_body_loads (@noinline) now returns panel_loads plus the moment about the reference point, for calc_forces! and for the body_aerodynamics results loop, which carried the same moment line. The POLAR_MATRICES testset also checks that the Jacobian does not depend on the ForwardDiff chunk size. Refs #360, #379 Co-Authored-By: Claude Opus 5 --- docs/src/private_functions.md | 1 + src/body_aerodynamics.jl | 9 ++++----- src/panel_aerodynamics.jl | 13 +++++++++++++ src/solver.jl | 9 ++++----- test/solver/test_forwarddiff.jl | 8 ++++++++ 5 files changed, 30 insertions(+), 10 deletions(-) diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index b6e76f60..68da6a8e 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -63,6 +63,7 @@ panel_force_directions panel_moment panel_couple_force panel_loads +panel_body_loads ``` ### Induced velocities diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 0d628876..b60a06a4 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -945,9 +945,10 @@ function calculate_results( cd_dist[i] += viscous.delta_cd c_span = viscous.c_span end - loads = panel_loads(axes, dirs, + arm = SVector{3}(panel.aero_center) - SVector{3}(reference_point) + loads = panel_body_loads(axes, dirs, dynamic_pressure(density, density, v_rel_dist[i]), - cl_dist[i], cd_dist[i], cm_dist[i]; c_span) + cl_dist[i], cd_dist[i], cm_dist[i], c_span, arm) force = loads.force va_panel = va_dist[i] @@ -980,11 +981,9 @@ function calculate_results( cd_prescribed_va[i] = drag_prescribed_va * inv_q_area cs_prescribed_va[i] = side_prescribed_va * inv_q_area - arm = SVector{3}(panel.aero_center) - SVector{3}(reference_point) - moment = loads.pitching_moment .* axes.y_airf .+ cross(arm, force) @inbounds for k in 1:3 f_body_3D[k, i] = force[k] - m_body_3D[k, i] = moment[k] + m_body_3D[k, i] = loads.body_moment[k] end end diff --git a/src/panel_aerodynamics.jl b/src/panel_aerodynamics.jl index 9cd2cd00..53b82ac1 100644 --- a/src/panel_aerodynamics.jl +++ b/src/panel_aerodynamics.jl @@ -198,3 +198,16 @@ along `y_airf`, as [`spanwise_flow_drag`](@ref) gives it. (c_span * q_dyn * chord) .* y_airf) return (; lift, drag, moment, force, pitching_moment=scale * width * moment) end + +""" + panel_body_loads(axes, dirs, q_dyn, cl, cd, cm, c_span, arm) + +[`panel_loads`](@ref) plus `body_moment`, the panel's moment about a reference point, +with `arm` the vector from that point to the aerodynamic center. +""" +@noinline function panel_body_loads(axes, dirs, q_dyn, cl, cd, cm, c_span, arm) + # Not inlined: Julia 1.12 miscompiles the inlined form for 10-partial Duals on Zen 4/5. + loads = panel_loads(axes, dirs, q_dyn, cl, cd, cm; c_span) + return (; loads..., + body_moment=loads.pitching_moment .* axes.y_airf .+ cross(arm, loads.force)) +end diff --git a/src/solver.jl b/src/solver.jl index 69dfb1b1..cd3143b5 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -452,19 +452,18 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; cd_dist[i] += viscous.delta_cd c_span = viscous.c_span end - loads = panel_loads(axes, dirs, + arm_vec = SVector{3, T}(panel.aero_center) - SVector{3, T}(reference_point) + loads = panel_body_loads(axes, dirs, dynamic_pressure(density, density, v_rel_dist[i]), - cl_dist[i], cd_dist[i], cm_dist[i]; c_span) + cl_dist[i], cd_dist[i], cm_dist[i], c_span, arm_vec) lift[i] = loads.lift drag[i] = loads.drag panel_moment_dist[i] = loads.moment force = loads.force - arm_vec = SVector{3, T}(panel.aero_center) - SVector{3, T}(reference_point) - moment = loads.pitching_moment .* axes.y_airf .+ cross(arm_vec, force) @inbounds for k in 1:3 solver.sol.f_body_3D[k, i] = force[k] - solver.sol.m_body_3D[k, i] = moment[k] + solver.sol.m_body_3D[k, i] = loads.body_moment[k] end arm = (moment_frac - 0.25) * panel.chord diff --git a/test/solver/test_forwarddiff.jl b/test/solver/test_forwarddiff.jl index 6fc3d4f6..d3c6aa59 100644 --- a/test/solver/test_forwarddiff.jl +++ b/test/solver/test_forwarddiff.jl @@ -104,5 +104,13 @@ end @info "POLAR_MATRICES jacobian norms" norm_fwd=norm(jac_fwd) norm_fd=norm(jac_fd) @test relative_error(jac_fd, jac_fwd) < 1e-4 + + @testset "the Jacobian does not depend on the ForwardDiff chunk size" begin + jac_chunk5, _, _ = VortexStepMethod.linearize( + matrix_solver, matrix_body, y_op; + theta_idxs=1:4, va_vec_idxs=5:7, omega_idxs=8:10, + aero_coeffs=true, backend=AutoForwardDiff(chunksize=5)) + @test relative_error(jac_chunk5, jac_fwd) < 1e-12 + end end end From 961455441578754f0b3a0bface97428d6b378c15 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Tue, 22 Sep 2026 14:50:50 +0200 Subject: [PATCH 4/6] Changelog: linearize's ForwardDiff Jacobian on Julia 1.12 and Zen 4/5 Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 22410bed..b826d7e7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -66,6 +66,9 @@ ### Fixed +- `linearize` with `AutoForwardDiff` gives the right Jacobian under Julia 1.12 on AMD + Zen 4/5 CPUs. With ten inputs, ForwardDiff's full chunk, Julia 1.12 miscompiled the + per-panel force assembly there and the Jacobian came out about 5% off. - `ELLIPTIC` initial circulation works on a body with more than one wing, where it threw an `ArgumentError`: each wing gets an ellipse over its own span, along its own `spanwise_direction` and centred on its own mid-span, also for a single wing off y = 0. From b8a8e4aeaadbac388bfa1992b2071ac18ff40b64 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Tue, 22 Sep 2026 15:56:32 +0200 Subject: [PATCH 5/6] Say which runs Julia 1.12 miscompiles the force assembly in: --check-bounds=yes Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 7 ++++--- src/panel_aerodynamics.jl | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index af38608a..9a1cdc99 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -77,9 +77,10 @@ ### Fixed -- `linearize` with `AutoForwardDiff` gives the right Jacobian under Julia 1.12 on AMD - Zen 4/5 CPUs. With ten inputs, ForwardDiff's full chunk, Julia 1.12 miscompiled the - per-panel force assembly there and the Jacobian came out about 5% off. +- `linearize` with `AutoForwardDiff` gives the right Jacobian under Julia 1.12 with + `--check-bounds=yes`, as `Pkg.test` runs, on some AMD Zen 4/5 CPUs. With ten inputs, + ForwardDiff's full chunk, Julia miscompiled the per-panel force assembly there and the + Jacobian came out about 5% off. - On a body whose wings span different directions, such as a wing and a vertical fin, `solve!`, `solve` and `linearize` take each panel's lift, drag and side directions from its own wing's `spanwise_direction`, not the first wing's. `solve` computes diff --git a/src/panel_aerodynamics.jl b/src/panel_aerodynamics.jl index 242c7df4..8f0063a3 100644 --- a/src/panel_aerodynamics.jl +++ b/src/panel_aerodynamics.jl @@ -204,7 +204,7 @@ end with `arm` the vector from that point to the aerodynamic center. """ @noinline function panel_body_loads(axes, dirs, q_dyn, cl, cd, cm, c_span, arm) - # Not inlined: Julia 1.12 miscompiles the inlined form for 10-partial Duals on Zen 4/5. + # Not inlined: Julia 1.12 --check-bounds=yes miscompiles it for Dual{10} on Zen 4/5. loads = panel_loads(axes, dirs, q_dyn, cl, cd, cm; c_span) return (; loads..., body_moment=loads.pitching_moment .* axes.y_airf .+ cross(arm, loads.force)) From 25768667697e2c14eaa8c35badadd858311ef8b3 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Tue, 22 Sep 2026 16:28:54 +0200 Subject: [PATCH 6/6] Name the LLVM 18 SLP-vectorizer bug behind the force-assembly miscompile Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 6 +++--- src/panel_aerodynamics.jl | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9a1cdc99..97aadbbc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -78,9 +78,9 @@ ### Fixed - `linearize` with `AutoForwardDiff` gives the right Jacobian under Julia 1.12 with - `--check-bounds=yes`, as `Pkg.test` runs, on some AMD Zen 4/5 CPUs. With ten inputs, - ForwardDiff's full chunk, Julia miscompiled the per-panel force assembly there and the - Jacobian came out about 5% off. + `--check-bounds=yes`, as `Pkg.test` runs, on AMD Zen 4/5 CPUs with AVX-512. With ten + inputs, ForwardDiff's full chunk, LLVM 18's SLP vectorizer miscompiled the per-panel + force assembly there (JuliaLang/julia#62368) and the Jacobian came out about 5% off. - On a body whose wings span different directions, such as a wing and a vertical fin, `solve!`, `solve` and `linearize` take each panel's lift, drag and side directions from its own wing's `spanwise_direction`, not the first wing's. `solve` computes diff --git a/src/panel_aerodynamics.jl b/src/panel_aerodynamics.jl index 8f0063a3..e2e0c67f 100644 --- a/src/panel_aerodynamics.jl +++ b/src/panel_aerodynamics.jl @@ -204,7 +204,7 @@ end with `arm` the vector from that point to the aerodynamic center. """ @noinline function panel_body_loads(axes, dirs, q_dyn, cl, cd, cm, c_span, arm) - # Not inlined: Julia 1.12 --check-bounds=yes miscompiles it for Dual{10} on Zen 4/5. + # Not inlined: Julia 1.12's LLVM SLP-vectorizes it wrongly for Dual{10} on AVX-512 Zen. loads = panel_loads(axes, dirs, q_dyn, cl, cd, cm; c_span) return (; loads..., body_moment=loads.pitching_moment .* axes.y_airf .+ cross(arm, loads.force))