Name linearize's locals and share the two-wing test setup

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Author: 1-Bort-1

src/body_aerodynamics.jl

@@ -260,13 +260,13 @@ which run over the unrefined sections of all wings in order; `nothing` leaves th
 unchanged. Call [`reinit!`](@ref) afterwards to update the panels.
 """
 function unrefined_deform!(body_aero::BodyAerodynamics, theta_angles, delta_angles)
-    first_section = 1
+    section_offset = 0
     for wing in body_aero.wings
-        wing_sections = first_section:first_section + wing.n_unrefined_sections - 1
+        section_idxs = section_offset .+ (1:wing.n_unrefined_sections)
         unrefined_deform!(wing,
-            isnothing(theta_angles) ? nothing : view(theta_angles, wing_sections),
-            isnothing(delta_angles) ? nothing : view(delta_angles, wing_sections))
-        first_section += wing.n_unrefined_sections
+            isnothing(theta_angles) ? nothing : view(theta_angles, section_idxs),
+            isnothing(delta_angles) ? nothing : view(delta_angles, section_idxs))
+        section_offset += wing.n_unrefined_sections
     end
     return nothing
 end

src/solver.jl

@@ -1177,16 +1177,19 @@ function _wing_with_eltype(wing::Wing{P, Float64}, ::Type{TD}) where {P, TD}
         wing.spanwise_distribution,
         PanelProperties{P, TD}(),
         MVector{3, TD}(wing.spanwise_direction),
-        Section{TD}[_section_with_eltype(s, TD) for s in wing.unrefined_sections],
-        Section{TD}[_section_with_eltype(s, TD) for s in wing.refined_sections],
+        Section{TD}[_section_with_eltype(section, TD)
+                    for section in wing.unrefined_sections],
+        Section{TD}[_section_with_eltype(section, TD)
+                    for section in wing.refined_sections],
         wing.remove_nan,
         wing.use_prior_polar,
         wing.billowing_percentage,
         TD(wing.crease_frac),
         copy(wing.refined_panel_mapping),
         copy(wing.refined_section_left_idx),
         Vector{TD}(wing.refined_section_weight),
-        Section{TD}[_section_with_eltype(s, TD) for s in wing.non_deformed_sections],
+        Section{TD}[_section_with_eltype(section, TD)
+                    for section in wing.non_deformed_sections],
         Vector{TD}(wing.theta_dist),
         Vector{TD}(wing.delta_dist),
         TD(wing.mass),
@@ -1261,7 +1264,7 @@ Returns `(jac, results, converged)` where `results` is `(F, M, moment_unrefined_
 or the corresponding coefficients when `aero_coeffs=true` — and `converged` is `false` (with a
 warning) if any internal solve missed the solver's tolerances.
 """
-function linearize(solver::Solver{P, U}, body_aero::BodyAerodynamics, y::Vector{T};
+function linearize(solver::Solver{<:Any, U}, body_aero::BodyAerodynamics, y::Vector{T};
         theta_idxs=1:4,
         delta_idxs=nothing,
         va_idxs=nothing,
@@ -1270,11 +1273,12 @@ function linearize(solver::Solver{P, U}, body_aero::BodyAerodynamics, y::Vector{
         backend = AutoForwardDiff(),
         fd_absstep::Float64=1e-8,
         fd_relstep::Float64=1e-8,
-        kwargs...) where {P, U, T}
+        kwargs...) where {U, T}
 
     for (name, idxs) in (("theta_idxs", theta_idxs), ("delta_idxs", delta_idxs))
         isnothing(idxs) || length(idxs) == U || throw(ArgumentError(
-            "Length of $name ($(length(idxs))) must match number of unrefined sections ($U)"))
+            "Length of $name ($(length(idxs))) must match number of unrefined sections " *
+            "($U)"))
     end
 
     n_failed = Ref(0)
@@ -1302,8 +1306,8 @@ function linearize(solver::Solver{P, U}, body_aero::BodyAerodynamics, y::Vector{
         end
 
         va = isnothing(va_idxs) ? MVector{3, TI}(body_aero_c._va) : y_in[va_idxs]
-        om = isnothing(omega_idxs) ? MVector{3, TI}(body_aero_c.omega) : y_in[omega_idxs]
-        set_va!(body_aero_c, va, om)
+        omega = isnothing(omega_idxs) ? MVector{3, TI}(body_aero_c.omega) : y_in[omega_idxs]
+        set_va!(body_aero_c, va, omega)
 
         solve!(solver_c, body_aero_c; kwargs...)
         solver_c.lr.converged || (n_failed[] += 1)
@@ -1319,13 +1323,13 @@ function linearize(solver::Solver{P, U}, body_aero::BodyAerodynamics, y::Vector{
         return nothing
     end
 
-    n_results = 3 + 3 + length(solver.sol.moment_unrefined_dist)
+    n_results = 3 + 3 + U
     jac = zeros(n_results, length(y))
     results = zeros(n_results)
-    be = backend === nothing ?
+    ad_backend = backend === nothing ?
         AutoFiniteDiff(absstep=fd_absstep, relstep=fd_relstep) : backend
-    prep = prepare_jacobian(calc_results!, results, be, y)
-    jacobian!(calc_results!, results, jac, prep, be, y)
+    prep = prepare_jacobian(calc_results!, results, ad_backend, y)
+    jacobian!(calc_results!, results, jac, prep, ad_backend, y)
     calc_results!(results, y)
     converged = n_failed[] == 0
     if !converged

test/body_aerodynamics/test_body_aerodynamics.jl

@@ -513,15 +513,23 @@ function solve_wings(wings)
 end
 
 """
-    linearize_wings(wings; kwargs...)
+    wing_pair(section_y, n_panels, offset)
 
-`linearize` at zero twist and deflection of the body [`solve_wings`](@ref) builds from
-`wings`, over the twist and the deflection of every unrefined section, then the inflow
-and the angular rate.
+Two `inviscid_wing`s of `n_panels` panels, the second shifted by `-offset` [m] in y.
 """
-function linearize_wings(wings; kwargs...)
-    n_sections = sum(wing -> wing.n_unrefined_sections, wings)
-    body_aero = BodyAerodynamics(wings; va=[10.0, 0.0, 1.0])
+function wing_pair(section_y, n_panels, offset)
+    return [inviscid_wing(section_y; n_panels),
+            inviscid_wing(section_y .- offset; n_panels)]
+end
+
+"""
+    linearize_body(body_aero; kwargs...)
+
+`linearize` of `body_aero` at zero twist and deflection over the twist and the deflection of
+every unrefined section, then the inflow and the angular rate.
+"""
+function linearize_body(body_aero; kwargs...)
+    n_sections = sum(wing -> wing.n_unrefined_sections, body_aero.wings)
     solver = Solver(body_aero; use_gamma_prev=false, rtol=1e-10)
     y0 = [zeros(2n_sections); body_aero.va; zeros(3)]
     return VortexStepMethod.linearize(solver, body_aero, y0;
@@ -538,9 +546,7 @@ end
     @test single.solver_status == FEASIBLE
 
     @testset "wings far apart each act as the isolated wing" begin
-        wings = [inviscid_wing(section_y; n_panels),
-                 inviscid_wing(section_y .- 1e4; n_panels)]
-        body_aero, sol = solve_wings(wings)
+        body_aero, sol = solve_wings(wing_pair(section_y, n_panels, 1e4))
 
         @test length(body_aero.panels) == 2n_panels
         @test sol.solver_status == FEASIBLE
@@ -552,9 +558,7 @@ end
     end
 
     @testset "wings one chord apart induce on each other" begin
-        wings = [inviscid_wing(section_y; n_panels),
-                 inviscid_wing(section_y .- (span + 1.0); n_panels)]
-        _, sol = solve_wings(wings)
+        _, sol = solve_wings(wing_pair(section_y, n_panels, span + 1.0))
         gamma = sol.gamma_distribution
 
         @test sol.solver_status == FEASIBLE
@@ -563,31 +567,30 @@ end
         @test sol.force[3] > 2single.force[3]
     end
 
-    n_sections = 2length(section_y)
+    n_wing_sections = length(section_y)
     @testset "linearize: theta of each wing moves that wing's sections" begin
-        wings = [inviscid_wing(section_y; n_panels),
-                 inviscid_wing(section_y .- 1e4; n_panels)]
-        jac, _, converged = linearize_wings(wings)
-        wing_rows = (7:6+length(section_y), 7+length(section_y):6+n_sections)
-        wing_columns = (1:length(section_y), length(section_y)+1:n_sections)
+        body_aero, _ = solve_wings(wing_pair(section_y, n_panels, 1e4))
+        jac, _, converged = linearize_body(body_aero)
+        first_sections = 1:n_wing_sections
+        second_sections = n_wing_sections+1:2n_wing_sections
+        own_first = jac[6 .+ first_sections, first_sections]
 
         @test converged
-        @test norm(jac[wing_rows[1], wing_columns[1]]) > 0
-        @test jac[wing_rows[2], wing_columns[2]] ≈ jac[wing_rows[1], wing_columns[1]] rtol=1e-4
-        @test norm(jac[wing_rows[1], wing_columns[2]]) < 1e-4norm(jac[wing_rows[1], wing_columns[1]])
-        @test norm(jac[wing_rows[2], wing_columns[1]]) < 1e-4norm(jac[wing_rows[1], wing_columns[1]])
+        @test norm(own_first) > 0
+        @test jac[6 .+ second_sections, second_sections] ≈ own_first rtol=1e-4
+        @test norm(jac[6 .+ first_sections, second_sections]) < 1e-4norm(own_first)
+        @test norm(jac[6 .+ second_sections, first_sections]) < 1e-4norm(own_first)
     end
 
     @testset "linearize: AutoForwardDiff matches AutoFiniteDiff" begin
-        wings = [inviscid_wing(section_y; n_panels),
-                 inviscid_wing(section_y .- (span + 1.0); n_panels)]
-        jac_fwd, _, fwd_converged = linearize_wings(wings)
-        jac_fd, _, fd_converged = linearize_wings(wings; backend=nothing,
+        body_aero, _ = solve_wings(wing_pair(section_y, n_panels, span + 1.0))
+        jac_fwd, _, fwd_converged = linearize_body(body_aero)
+        jac_fd, _, fd_converged = linearize_body(body_aero; backend=nothing,
             fd_absstep=1e-6, fd_relstep=1e-6)
 
         @test fwd_converged
         @test fd_converged
-        @test norm(jac_fwd[:, 1:n_sections]) > 0
+        @test norm(jac_fwd[:, 1:2n_wing_sections]) > 0
         @test jac_fwd ≈ jac_fd rtol=1e-4
     end
 end