Keep calc_forces!' per-panel loads out of line so Julia 1.12 gets linearize's ForwardDiff Jacobian right on AVX-512 Zen

CHANGELOG.md

@@ -8,6 +8,10 @@
   the right tip, z = x × y up. The separate `KB` definition, whose `Z = Y × X` pointed
   down, is gone; the solver's frame is unchanged. CL, CD and CS are documented as
   wind-axis coefficients, `cfx`, `cfy` and `cfz` as the body-axis ones.
+- `linearize` with `AutoForwardDiff` gives the right Jacobian under Julia 1.12 with
+  `--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.
 
 ## VortexStepMethod v6.0.0 2026-09-23
 

docs/src/private_functions.md

@@ -64,6 +64,7 @@ prescribed_va_directions
 panel_moment
 panel_couple_force
 panel_loads
+panel_body_loads
 ```
 
 ### Induced velocities

src/body_aerodynamics.jl

@@ -953,9 +953,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]
@@ -981,11 +982,9 @@ function calculate_results(
             cd_prescribed_va[i] = drag_prescribed_va * inv_q_area
             cs_prescribed_va[i] = dot(force, wing_dirs.dir_side) * 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
     end

src/panel_aerodynamics.jl

@@ -196,3 +196,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'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))
+end

src/solver.jl

@@ -449,19 +449,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

test/solver/test_forwarddiff.jl

@@ -97,5 +97,13 @@ relative_error(jac, reference) = maximum(abs.(jac .- reference)) / maximum(abs,
 
         @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(
+                ram_solver, ram_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