[ocayon]

Summary

Two changes, both behaviour-preserving (no physics change):

1. The quarter-chord force law becomes a public API and xp kernels, so CasADi callers (AWETrim's CasADi trim) can use VSM's formula instead of restating it. compute_results calls the same kernels for its own per-panel force, so the formula lives in one place.
2. Allocation-free numba filament kernels. The AIC assembly is about 70× faster and a casadi_newton solve about 1.8× faster.

New public functions (for the AWETrim follow-up)

BodyAerodynamics

compute_quarter_chord_AIC_matrices(core_radius_fraction, va_norm_array, va_unit_array)
    -> (AIC_x, AIC_y, AIC_z)            # each (n, n) = AIC_LLT + AIC_curved_correction

• This is the induction of the force-direction flow at the quarter chords: v_qc = va + AIC_qc @ gamma.
• It uses the existing per-shape cache of the curved bound-vortex correction, which is rotated under rigid motion.
• update_effective_angle_of_attack_if_VSM now uses it, so there is one source.

VSM.core.panel_kernels

Every function takes xp=numpy (the default) or xp=casadi. Conventions: (n, 3) rows; (n,) vectors, which are (n, 1) columns in CasADi.

relative_flow(v_rel, x_airf, y_airf, z_airf, xp=np) -> (alpha, umag)
quarter_chord_velocity(va, gamma, AIC_x, AIC_y, AIC_z, xp=np) -> v (n, 3)     # va + AIC gamma
direction_angle(v, x_airf, y_airf, xp=np) -> alpha_dir                        # atan2(v.x_airf, v.y_airf)
section_force_directions(alpha_dir, x_airf, y_airf, z_airf, xp=np) -> (dir_flow, dir_lift)
section_force_magnitudes(umag, cl, cd, cm, chord, rho) -> (lift, drag, moment)  # per unit span
section_force_vectors(alpha, umag, cl, cd, cm, chord, x_airf, y_airf, z_airf, rho,
                      alpha_dir=None, xp=np) -> (lift_2d, drag_2d, moment, dir_flow)
panel_forces(alpha, umag, cl, cd, cm, chord, width, x_airf, y_airf, z_airf, rho,
             alpha_dir=None, xp=np) -> (force (n, 3), moment_local (n, 3))
panel_moment_local(moment, width, z_airf, xp=np) -> (n, 3)
attached_trailed_forces(gamma, va_panels, v_induced, legs, rho, xp=np) -> (n+1, 3)
# row helpers: hstack_columns, cross_rows, dot_rows, norm_rows, scale_rows, unit_rows

• The force magnitudes come from the polar at the 3/4-chord alpha with q = 0.5 rho umag^2.
• The force directions come from alpha_dir:
  • alpha_dir=None uses alpha itself: the is_aoa_corrected=False law.
  • direction_angle(quarter_chord_velocity(va, gamma, *AIC_qc)) gives the True law (TAT3).
• The viscous drag correction stays NumPy-only. It is added on top in compute_results, as before.

Solver

casadi_panel_force_function() -> casadi.Function
  inputs : v_rel (n,3), v_qc (n,3), x_airf, y_airf, z_airf (n,3), chord, width (n), rho
  outputs: force (n,3), moment_local (n,3), alpha, alpha_dir, umag, cl, cd, cm (n)

• It is cached per polar set and uses the same polar lookups as the casadi_newton residual. Those lookups are factored into _casadi_polar_lookups; the residual function's inputs, outputs and graph are unchanged.
• It takes v_qc rather than the matrices, so the graph stays O(n). The caller forms the dense product, the same design as the Newton residual taking v_rel.
• It returns the panel forces without the attached-trailed force. That force is available as panel_kernels.attached_trailed_forces.

What AWETrim should call

1. Per wake pass, freeze numerically:
   • AIC_qc = body.compute_quarter_chord_AIC_matrices(core_radius_fraction, va_norm, va_unit);
   • the 3/4-chord AIC matrices it already freezes.
2. In the graph:
   • v_rel = panel_kernels.quarter_chord_velocity(va, gamma, *AIC_cp, xp=ca);
   • v_qc = panel_kernels.quarter_chord_velocity(va, gamma, *AIC_qc, xp=ca) for is_aoa_corrected=True, or v_qc = v_rel for False.
3. force, moment_local, ... = solver.casadi_panel_force_function()(v_rel, v_qc, x_airf, y_airf, z_airf, chord, width, rho).
4. For the attached-trailed force, use panel_kernels.attached_trailed_forces(...). Its xp conventions match AWETrim's panel_kernels.py.

The row helpers and relative_flow / panel_forces / attached_trailed_forces keep the names and conventions of AWETrim's panel_kernels.py. Two differences:

• panel_forces has no spanwise_direction argument (it was unused there) and takes an optional alpha_dir.
• panel_forces keeps VSM's operation order (cm 0.5 rho U^2 c^2) so it matches compute_results exactly.

No-change regression

Reference. tests/BodyAerodynamics/data/panel_kernel_reference.npz holds every numeric output of Solver.solve (576 arrays), written with VSM at 1201f6f, the original code before both commits. It covers six cases, each with is_aoa_corrected True and False:

• a straight wing, a 30° swept wing and a circular-arc wing;
• the LEI V3 kite three ways: plain; with the viscous drag correction and no attached-trailed force; and with a yaw rate and an offset reference point.

Largest deviation from 1201f6f, relative to the quantity's own scale, or to the case's force/moment scale for entries that are zero by symmetry:

case	True	False
straight	1.3e-15	2.4e-15
swept	7.4e-15	1.5e-14
arc	4.2e-15	2.4e-15
V3	7.6e-15	1.1e-14
V3, viscous correction, no AT	8.1e-15	1.1e-14
V3, yaw rate, offset reference point	6.1e-15	6.7e-15

Other checks:

• CasADi kernels (numeric DM inputs) equal the NumPy kernels to 1e-12.
• A CasADi graph with gamma, va and both sets of induction matrices as MX parameters reproduces the compute_results panel forces at a converged gamma: force 2e-16, local moment 8e-15, for both laws.

Note on the CasADi comparison. solve passes compute_results the alpha/Umag of the circulation loop's last iterate. Those lag the returned gamma by one iteration, so they agree with it only to allowed_error. The graph test therefore compares against compute_results evaluated at one consistent gamma. This is existing behaviour and is not changed here. A CasADi trim that evaluates the force at its own gamma will differ from solve() by O(allowed_error).

Filament kernels (utils.py)

Cause. The numba kernels carried vectors as length-3 arrays. Every temporary (np.cross, r1 - r0, np.zeros(3)) was a heap allocation, and those allocations were about 95 % of the cost.

Change. The kernels now carry 3-tuples, which numba keeps on the stack. The formulas, branches and summation order are unchanged. error_model="numpy" keeps NumPy's inf/nan semantics for zero denominators, where numba's default would raise ZeroDivisionError on scalars, so degenerate geometry behaves as before. The existing equivalence tests against Filament.py (rtol 1e-12) pass.

	before	after
per filament evaluation	1.1 µs	15-19 ns
AIC, 36 / 100 panels	7.3 / 57 ms	0.10 / 0.93 ms
V3 casadi_newton solve, 36 panels	44 ms	24 ms

This also speeds up assemble_curved_bound_correction_AIC and induced_velocity_at_points (the attached-trailed force).

Post-processing and Newton-call overhead (third commit)

• Trefftz drag: compute_trefftz_plane_induced_drag is a far-wake diagnostic computed on every solve. Its O(n²) Python loop is now vectorised: 1.3 → 0.2 ms at 36 panels.
• Projected area: Wing.compute_projected_area is vectorised: 1.3 → 0.1 ms. It also no longer normalises the caller's z_plane_vector in place.
• Newton residual calls: gamma_loop_casadi_newton converts the inputs that are constant within a solve to CasADi matrices once, and calls a cached six-output variant of the residual function (same graph): ~185 → ~77 µs per evaluation.

Outputs move at round-off only (≤ 1.5e-14 against 1201f6f). The test counts are unchanged.

V3, 36 panels, best of 7 sweeps (machine under other load), casadi_newton: 50.7 ms at 1201f6f → 23.8 ms after commit 2 → 17.0 ms after commit 3.

Tests

Full suite: 201 passed, 5 skipped. That is the 183 passed, 5 skipped at 1201f6f plus 18 new tests in tests/BodyAerodynamics/test_panel_kernels.py.

Written by Claude, reviewed by Oriol Cayon

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has insufficient remaining quota. Try requesting a review at a lower effort level.

An organization owner can choose the repository organization as the billing source for reviews requested by members with a paid Copilot license. Learn about Copilot code review billing. Review the Copilot code review billing settings.