linearize and make_dual_shadow take a BodyAerodynamics with several wings

CHANGELOG.md

@@ -2,6 +2,11 @@
 
 ## Unreleased
 
+### Added
+
+- `linearize` takes a `BodyAerodynamics` with more than one wing; `theta_idxs` and
+  `delta_idxs` then run over the unrefined sections of all wings in order.
+
 ### Changed
 
 - Requires Julia 1.12 or 1.13; 1.10 and 1.11 keep resolving v5.1.1.

src/body_aerodynamics.jl

@@ -252,6 +252,25 @@ function calculate_stall_angle_list!(stall_angles::AbstractVector,
     return nothing
 end
 
+"""
+    unrefined_deform!(body_aero::BodyAerodynamics, theta_angles, delta_angles)
+
+Deform each wing of `body_aero` by its entries of `theta_angles` and `delta_angles` [rad],
+which run over the unrefined sections of all wings in order; `nothing` leaves that angle
+unchanged. Call [`reinit!`](@ref) afterwards to update the panels.
+"""
+function unrefined_deform!(body_aero::BodyAerodynamics, theta_angles, delta_angles)
+    section_offset = 0
+    for wing in body_aero.wings
+        section_idxs = section_offset .+ (1:wing.n_unrefined_sections)
+        unrefined_deform!(wing,
+            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
+
 """
     reinit!(body_aero::BodyAerodynamics; init_aero, va, omega, refine_mesh, recompute_mapping, sort_sections)
 

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),
@@ -1213,10 +1216,8 @@ buffers are freshly allocated as `TD`-typed.
 function make_dual_shadow(solver::Solver{P, U, Float64},
                           body_aero::BodyAerodynamics{P, W, Float64},
                           ::Type{TD}) where {P, U, W, TD}
-    length(body_aero.wings) == 1 || throw(ArgumentError(
-        "make_dual_shadow currently supports body_aero with one wing"))
-    wing_d = _wing_with_eltype(body_aero.wings[1], TD)
-    body_aero_d = BodyAerodynamics([wing_d];
+    wings_d = [_wing_with_eltype(wing, TD) for wing in body_aero.wings]
+    body_aero_d = BodyAerodynamics(wings_d;
         va = MVector{3, TD}(body_aero._va),
         omega = MVector{3, TD}(body_aero.omega),
     )
@@ -1251,8 +1252,8 @@ end
               backend=AutoForwardDiff(), kwargs...)
 
 Jacobian of aerodynamic outputs w.r.t. control and kinematic inputs at `y`. Each `*_idxs`
-selects which entries of `y` map to twist angles (one per unrefined section), trailing-edge
-deflections (one per unrefined section), apparent wind `(vx, vy, vz)`, and angular rate
+selects which entries of `y` map to twist angles and trailing-edge deflections (one per
+unrefined section, over all wings in order), apparent wind `(vx, vy, vz)`, and angular rate
 `(ωx, ωy, ωz)` respectively.
 
 `backend` accepts any `DifferentiationInterface` backend; `AutoForwardDiff()` (the default)
@@ -1263,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, 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,
@@ -1272,22 +1273,12 @@ function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T};
         backend = AutoForwardDiff(),
         fd_absstep::Float64=1e-8,
         fd_relstep::Float64=1e-8,
-        kwargs...) where T
+        kwargs...) where {U, T}
 
-    !(length(body_aero.wings) == 1) && throw(ArgumentError("Linearization only works for a body_aero with one wing"))
-    wing = body_aero.wings[1]
-
-    # Validate that theta_idxs and delta_idxs match the number of unrefined sections
-    if !isnothing(theta_idxs) && wing.n_unrefined_sections > 0
-        length(theta_idxs) != wing.n_unrefined_sections && throw(ArgumentError(
-            "Length of theta_idxs ($(length(theta_idxs))) must match number of unrefined sections ($(wing.n_unrefined_sections))"))
-    end
-    if !isnothing(delta_idxs) && wing.n_unrefined_sections > 0
-        length(delta_idxs) != wing.n_unrefined_sections && throw(ArgumentError(
-            "Length of delta_idxs ($(length(delta_idxs))) must match number of unrefined sections ($(wing.n_unrefined_sections))"))
-    end
-    if wing.n_unrefined_sections == 0 && (!isnothing(theta_idxs) || !isnothing(delta_idxs))
-        throw(ArgumentError("Cannot use theta_idxs or delta_idxs when wing has no unrefined sections"))
+    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)"))
     end
 
     n_failed = Ref(0)
@@ -1298,27 +1289,25 @@ function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T};
         if TI === Float64
             body_aero_c = body_aero
             solver_c = solver
-            wing_c = wing
         else
             shadow = shadow_ref[]
             if shadow === nothing || eltype(shadow[1]._va) !== TI
                 shadow_ref[] = make_dual_shadow(solver, body_aero, TI)
             end
             body_aero_c, solver_c = shadow_ref[]
-            wing_c = body_aero_c.wings[1]
         end
 
         @views theta_angles = isnothing(theta_idxs) ? nothing : y_in[theta_idxs]
         @views delta_angles = isnothing(delta_idxs) ? nothing : y_in[delta_idxs]
 
         if !isnothing(theta_angles) || !isnothing(delta_angles)
-            VortexStepMethod.unrefined_deform!(wing_c, theta_angles, delta_angles; smooth=false)
-            VortexStepMethod.reinit!(body_aero_c; init_aero=false)
+            unrefined_deform!(body_aero_c, theta_angles, delta_angles)
+            reinit!(body_aero_c; init_aero=false)
         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)
@@ -1334,13 +1323,13 @@ function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T};
         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

@@ -512,6 +512,32 @@ function solve_wings(wings)
     return body_aero, solve!(Solver(body_aero), body_aero)
 end
 
+"""
+    wing_pair(section_y, n_panels, offset)
+
+Two `inviscid_wing`s of `n_panels` panels, the second shifted by `-offset` [m] in y.
+"""
+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;
+        theta_idxs=1:n_sections, delta_idxs=n_sections+1:2n_sections,
+        va_idxs=2n_sections+1:2n_sections+3, omega_idxs=2n_sections+4:2n_sections+6,
+        kwargs...)
+end
+
 @testset "solve! on a two-wing body" begin
     n_panels = 6
     section_y = [2.0, 0.0, -2.0]
@@ -520,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
@@ -534,14 +558,39 @@ 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
         @test gamma ≈ reverse(gamma) rtol=1e-6
         @test gamma[n_panels] > 1.05single.gamma_distribution[n_panels]
         @test sol.force[3] > 2single.force[3]
     end
+
+    n_wing_sections = length(section_y)
+    @testset "linearize: theta of each wing moves that wing's sections" begin
+        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(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
+        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:2n_wing_sections]) > 0
+        @test jac_fwd ≈ jac_fd rtol=1e-4
+    end
 end