Skip to content

The default solve now gives the same results as the Python VSM v2.3.0 - #415

Draft
1-Bort-1 wants to merge 9 commits into
agent/94-remove-results-dictfrom
agent/410-vortexstepmethod-jl-has-drifted-from-the
Draft

1-Bort-1 wants to merge 9 commits into
agent/94-remove-results-dictfrom
agent/410-vortexstepmethod-jl-has-drifted-from-the

Conversation

@1-Bort-1

@1-Bort-1 1-Bort-1 commented Nov 5, 2026 •

Copy link
Copy Markdown
Member

Julia was translated from Python v1.0.0, and Python's v2.3.0 changed its default results. On swept and curved wings, which is every kite, the two had drifted several percent apart. With this branch, Julia matches Python's circulation, angles, per-panel forces and moments, totals and coefficients to six digits on nine reference cases generated by running the Python code, with both solved to a tight stop. At the default stop, circulation, angles and coefficients stay within 1e-5 of those results.

V3 kite at 8 degrees, spanwise circulation: Julia before and after against Python v2.3.0, same geometry and polars

wing, 8° CL before CL now = Python CD before CD now = Python
straight 0.6906 0.6910 0.0186 0.0151
swept, tapered 0.8456 0.8256 0.0364 0.0299
arc 0.5947 0.5987 0.0354 0.0325
V3 kite 1.0509 1.0992 0.4910 0.4575

What moved: circulation now follows the effective speed. The chord and chordwise axis are measured square to the bound vortex. Forces turn with the flow at the quarter chord, including a correction for a curved bound vortex, and drag runs along that flow. The force on the chordwise legs is on by default and shared per section boundary. Each wake follows its own panel's inflow, and reference pressure and directions leave out body rates. The defaults also change to 5000 iterations and a relaxation factor of 0.8 of the stability limit, so the curved reference wing that used to diverge now converges.

The same default change ends the lost induced drag of #318: the straight elliptic wing of aspect ratio 4 from Gaunaa, Li & Pirrung 2026 now reaches an Oswald efficiency of 0.9947, against 0.7295 with the quarter-chord force direction off.

Agreeing with Python does not make swept wings right. Against that paper's vortex-lattice results, Oswald efficiency is off by −8.2 % at 15 % sweep, −8.0 % at 10 % sweep with 25 % anhedral and −11.1 % at 25 % sweep, where the solver does not converge. The paper's own model is off by +5.6 %, −0.5 % and +1.1 %. Purely curved wings are off by ±2.4 %. These were measured on the previous head by the task behind #400, whose card plots them, and the chord direction on the swept wings is an assumption. #416 tracks the gap; closing it may mean leaving Python's rule, which is the owner's choice.

Anyone relying on the defaults gets new numbers, including KiteModels wherever its settings do not name these options. Setting the two force options to false and naming a relaxation factor restores those switches, but not the geometry or circulation changes.

Decisions I made:

  • Convergence is still measured on the unrelaxed residual at 1e-5, because it does not depend on the relaxation factor. Python's relaxed step at 1e-6 is similar in size on most wings and looser on the V3 kite. Julia's stop leaves circulation within 5.5e-6 and cl and cd within 5.4e-6 of the converged result.
  • Warm starting from the previous circulation stays on, because simulation loops rely on it. Past stall this can pick a different solution from Python, which starts from zero.
  • Sideslip from angles keeps Julia's unit-length apparent wind. Python's is not unit length (1.00045 at 10° and 10°), so the reference cases pass the inflow as a vector.
  • With several wings the curved correction acts only within each wing; Python refuses more than one wing.
  • The correction is kept per shape and rebuilt after any shape change. Python's reuse under rigid rotation is not ported.
  • The correction is not rebuilt when the inflow speed changes, because that would cost a rebuild per simulation step. Speed only sets the size of the trailing core, and rebuilds at 1, 10 and 100 m/s are identical to the last bit on all nine cases. A point that did fall inside a core would keep the speed of the last rebuild.
  • Drag beyond a polar table is now held, as in Python, for vector polars only. Matrix polars have no Python counterpart and are unchanged.
  • Inside a trailing vortex's core the speed is no longer limited, a deliberate match to Python: closer than 7 micrometres to the line at 10 m/s, the speed is the value on the core boundary, 1.85e5 m/s at 1e-7 m for unit circulation. The nine reference cases still match.

Found on the way: faster convergence exposed that unrefined section coefficients were plain means, so they are now area-weighted, with the same code as #402. In that code the section chord becomes area over width, and the section moment becomes the sum of its panels rather than their mean. That moment, and the linearize outputs and Jacobian rows taken from it, now grow with the panels per section. A semi-infinite filament now gives zero rather than NaN on its own axis, as #408 also does. It also counts a point within rounding of its line as on the line, by the rule the finite filament uses; such a point got 1.85e9 m/s at 1e-15 m.

This covers default results only, so the issue stays open. The missing features are queued as parts of it. Left out: bridle drag, which stays in KiteModels as ruled on #113, and the other items declined there (flat area, rotate, cosine circulation, configurable aerodynamic-centre fractions) unless the owner says otherwise.

To act on: stacked on #392, which merges first. It takes a different route from #400, whose rule is not Python's, and from #408's core scaling. It conflicts with #402 in the solver, where the aggregation code is identical. SymbolicAWEModels calls the drag-direction helper with an argument that is now gone, so it needs a one-line change when it raises its bound to the release #392 forces.

Verification

  • Reproduced first, unchanged code against Python: V3 CL 1.0509 vs 1.0992, CD 0.4910 vs 0.4575, chord off by up to 0.33 m
  • test/verification/test_python_reference.jl: green on nine cases (straight, swept, arc, sideslip, LLT, viscous, body rates, distributed inflow, V3). Red before: run through a script on the unchanged code, since the test itself names functions this branch adds
  • test/verification/test_verification.jl: the aspect-ratio-4 test is green, 24/24 in the file; red with the quarter-chord force direction disabled (0.7295, 23/24)
  • Green: 20 affected test files, including test/bench.jl, test_verification.jl, test_curved_bound_correction.jl, test_forwarddiff.jl; the six a later filament change touches (test_bound_filament.jl, test_panel.jl, test_curved_bound_correction.jl, test_python_reference.jl, test_forwarddiff.jl, test/bench.jl) rerun on the previous head; on this head the ten files the last fixes touch, and the docs
  • Docs build clean in a fresh session · repo has no REUSE lint
  • Local full suite: running on this head (Julia 1.13.3, one cell of the matrix) · GitHub CI: this draft runs the Documentation job only, green on this head; the Julia matrix starts once it leaves draft
  • test_python_reference.jl: the default stop is within 1e-5 on all nine cases. test_semi_infinite_filament.jl: rounding off the line gives zero, red before at three distances. test_forwarddiff.jl: the stability limit is type-stable on dual numbers, red before (Union{Float64, Dual})
  • test_curved_bound_correction.jl: rebuilding the correction allocates under 2 kB, red before the fix, when a V3 solve after a shape change allocated 1,536,208 B at 18 panels
  • Allocations, V3 at 18 panels: a solve after a shape change 1744 B, against 208 B with the two force options off; a repeat solve 208 B either way; calc_forces! 0 B
  • Iterations from a cold start, new defaults against the two force options off at the old relaxation factor 0.03: straight 46 vs 322, swept 55 vs 294, arc 39 vs 300, V3 315 vs 223
  • Benchmark: not run. The box stayed above its load limit for the full 8-hour wait (load 34 of 20 cores at the end), so no timing is claimed here
  • Tests that follow the change, same system: thesis geometry expectations take the square chord and frame, same tolerance (Python loosened its own to 1e-3). Pinned rectangular-wing values move 0.1 % (force x −117.965 → −117.838) with the circulation map. Tests that compare against section-only formulas now set the leg force off. cl is checked against the control-point angle. The NaN assertion on the filament axis is now zero. One trim call tightens its angle tolerance to 1e-6, since its 1e-5 assertion sat at 1.17e-5.
  • Reference: Python a2a91cb, numba 0.68.1, numpy 2.5.4; regenerate with test/data/python_reference/generate_reference.py
  • Risk: simulation loops that change the shape every step now rebuild the correction every step, 148 filament evaluations per panel pair (47,952 at 18 panels). A solve without a named relaxation factor also makes 91 lift lookups per panel to find the limit (1,638 at 18 panels, no allocation). Neither time is measured.

Closes #318

Opened by 1-Bort-1, an AI agent working for @1-Bart-1.
Refs #410 · task VortexStepMethod.jl-410

1-Bort-1 and others added 4 commits November 5, 2026 00:54
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Once a distance below 1e-16 has returned, the core-axis test needs a
distance below 1e-24 * 4 alpha nu / va, so an inflow under 7e-13 m/s. A
solve refuses zero inflow before it builds the influence matrices.

The docstring named the wrong function and arguments; it now gives the
right ones and says what happens inside the core.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
From the #318 branch: Gaunaa et al. 2026's straight elliptic wing reaches
1 with the quarter-chord flow turning the force, and stays below 0.8
without it.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@1-Bort-1

1-Bort-1 commented Nov 5, 2026 •

Copy link
Copy Markdown
Member Author

Local full suite: PASS (19 min, Julia 1.13.3, one cell of the matrix)

@1-Bort-1 1-Bort-1 added agent:queued Agent task state agent:running Agent task state and removed agent:ci Agent task state agent:queued Agent task state labels Nov 5, 2026
1-Bort-1 and others added 2 commits November 5, 2026 09:36
reinit! of a bound filament allocated a scratch vector on every call that did not
pass one. The curved-bound correction reinitialises its scratch filament 148 times per
panel pair, so a default solve after a shape change allocated 1.5 MB at 18 panels. The
filament's own r0 holds the same difference, so the argument goes.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@1-Bort-1 1-Bort-1 added agent:ci Agent task state agent:queued Agent task state agent:running Agent task state agent:review Agent task state and removed agent:running Agent task state agent:ci Agent task state agent:queued Agent task state labels Nov 5, 2026

@1-Bort-1 1-Bort-1 left a comment

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Independent review (advisory)

Verdict: REQUEST CHANGES · 5 inline, 0 off the diff

Good

  • Default change follows #410's agreement goal: test_python_reference.jl compares nine Python cases, with rtol 1e-6 on circulation and angles and 1e-6 of scale on forces and totals (read, not run).
  • Influence matrices are pinned to Python at rtol 1e-10 on the swept and arc cases, covering the wake and curved-correction changes.
  • Scratch vectors don't collide: the filament kernels index work_vectors 1, 3, 4 and 5, while calculate_AIC_matrices! holds 8-11 (grep only).
  • attached_trailed_loads! nets circulation per boundary and gives each tip its whole load, consistent with its docstring (read).
  • Old paths are deleted rather than kept: panel_force_moment, the va_magw/v_acrossz buffers, the shared wake vector and the spanwise argument of panel_force_directions.
  • The card states the known swept-wing gap (#416), the downstream break in SymbolicAWEModels, and that the benchmark was not run.

Not good

  • src/solver.jl:50 — moment_unrefined_dist changes from a per-section mean to a sum, but the card and CHANGELOG line 34 name only coefficients and angles. linearize returns this field (solver.jl:1339), so its moment outputs and Jacobian rows silently scale with panels per section for downstream users, breaking the §5 changelog rule.
  • src/solver.jl:200 — The default rtol stays 1e-5 on the unrelaxed residual, departing from #410's 'the two versions should agree'. The reference cases run both codes at 1e-10 (generate_reference.py:101), so the claim that default results agree to six digits is not shown at the default stop.
  • src/filament.jl:226 — This on-axis guard is an absolute 1e-16, below coordinate rounding at metre scale, while the finite kernel at line 89 uses a relative 1e-10. A point on a wake line with about 1e-15 rounding gets a core-boundary speed estimated near 1e9 m/s here but zero there, so one rule exists in two forms (§2).
  • src/solver.jl:731 — something evaluates its arguments eagerly, so relaxation_factor_limit runs 91 cl lookups per panel on every solve_base!, even with a named factor. Simulation loops pay this each step, including on the card's restore path, and the card's risk note omits it.
  • src/body_aerodynamics.jl:651 — The cache key is shape and core radius, but the correction is built from va_dist (line 631), which sets the trailing core size, and set_va! does not invalidate it. A solve can reuse a correction from an earlier speed, so results can depend on history whenever a point falls inside a core, though that core is millimetre-scale and the effect is likely small.
  • Not run: my Julia measurements were refused here, so the costs and effects below come from reading the code.
  • Warm start and unit-length sideslip remain #410 departures from Python; the sideslip reference case passes the inflow as a vector, so sideslip from angles is not compared.
  • CHANGELOG lists relaxation_factor_limit under Added as if public, but it is not exported and is documented in private_functions.md.
  • chord_unrefined_dist accumulates area under a chord name before being converted at solver.jl:593, and that chord change from plain mean to area/width is not named either.
  • test_body_aerodynamics.jl:235 is 93 characters against the 92 limit, and lines 228-229 are a two-line comment block where §3 allows one line.
  • relaxation_factor_limit seeds ratio and steepest as Float64, so the Dual shadow path may be type-unstable (not measured).
  • reference.yaml adds 4250 generated lines; this is tolerable as a fixture with its generator, but its keep_matrices rows are the bulk of later regeneration diffs.

claude, rubric CLEAN_CODE.md. A different lab from the implementer
on purpose: a reviewer sharing its blind spots would not flag its mistakes.

Comment thread src/solver.jl
- `cl_unrefined_dist`::MVector{U, Float64}: Averaged lift coefficients for unrefined sections [-]
- `cd_unrefined_dist`::MVector{U, Float64}: Averaged drag coefficients for unrefined sections [-]
- `cm_unrefined_dist`::MVector{U, Float64}: Averaged airfoil moment coefficients for unrefined sections [-]
- `moment_unrefined_dist`::MVector{U, Float64}: Pitching moment per unrefined section,

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

MAJOR: moment_unrefined_dist changes from a per-section mean to a sum, but the card and CHANGELOG line 34 name only coefficients and angles. linearize returns this field (solver.jl:1339), so its moment outputs and Jacobian rows silently scale with panels per section for downstream users, breaking the §5 changelog rule.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Named in 05ae687: CHANGELOG and card now say the section moment is a sum, and that linearize's moment outputs and Jacobian rows grow with panels per section; the sum stays, since moment is additive and test_unrefined_dist pins it.

Comment thread src/solver.jl Resolved
Comment thread src/filament.jl Resolved
Comment thread src/solver.jl Resolved
Comment thread src/body_aerodynamics.jl Resolved
@1-Bort-1 1-Bort-1 added agent:queued Agent task state agent:running Agent task state and removed agent:review Agent task state agent:queued Agent task state labels Nov 5, 2026
A point within coordinate rounding of a semi-infinite wake line got a core
boundary speed of 1.85e9 m/s from an absolute 1e-16 guard, where the finite
kernel gave zero from a relative 1e-10 one. Both now share the relative rule.

Skip the stability limit's lookups when a relaxation factor is named, type its
seeds so the Dual path is stable, and pin the default stop to within 1e-5 of the
converged Python result. Name the moment sum and the chord change in the
changelog, and take the private limit out of Added.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@1-Bort-1 1-Bort-1 added agent:queued Agent task state agent:running Agent task state agent:ci Agent task state and removed agent:running Agent task state agent:queued Agent task state labels Nov 5, 2026

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

agent:ci Agent task state

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant