Check POLAR_MATRICES ForwardDiff against finite differences on affine tables, not NeuralFoil's ram-air polars

CHANGELOG.md

@@ -77,6 +77,10 @@
 
 ### Fixed
 
+- `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.
 - 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

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

@@ -5,6 +5,24 @@ using Test
 
 relative_error(jac, reference) = maximum(abs.(jac .- reference)) / maximum(abs, reference)
 
+# 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(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)]
+        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 +75,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()
+        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,18 +91,26 @@ 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
 
         @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