bmad-sim / bmad-sim/BeamTracking.jl
SciBmadStandard integrates combined-function bends in a single step, giving silently wrong tunes
- Dominant language
- Julia
- Stars
- 8
- Forks
- 8
- Avg merge
- 16h 2m
- Merged PRs (30d)
- 6
Description
`SciBmadStandard` falls back to `Symplectic()` — i.e. `n_steps = 1` — whenever no closed-form map is available, with no dependence on how strong the element is. For a strong combined-function bend this produces a map that is wrong by a large factor but still exactly symplectic, so nothing downstream flags it and the user gets a plausible-looking, badly wrong answer.
## Reproducer
A CESR dipole-quadrupole: `L = 2.3475 m`, `g = 0.0318633 m⁻¹`, `Kn1 = -0.438171 m⁻²`, i.e. `√|k₁|·L = 1.554 rad` of betatron phase inside one element.
```julia
using GTPSA, Beamlines, BeamTracking, LinearAlgebra
const L, G, K1 = 2.3475, 0.0318633057911049677, -0.438170844695771
const d = Descriptor(6, 1)
function themap(tm)
ele = SBend(L = L, g_ref = G, Kn0 = G, Kn1 = K1)
isnothing(tm) || (ele.tracking_method = tm)
bl = Beamline([ele]; pc_ref = 6e9, species_ref = Species("positron"))
b = Bunch(collect(transpose(@vars(d))))
track!(b, bl)
return GTPSA.jacobian(b.coords.v)
end
# analytic combined-function sector bend
kx, ky = K1 + G^2, -K1
sx, sy = sqrt(-kx), sqrt(ky)
println("analytic: kx = ", -sx*sinh(sx*L), " ky = ", sy*sin(sy*L))
M = themap(nothing)
println("default: kx = ", -M[2,1], " ky = ", -M[4,3])
for n in (1, 2, 4, 8, 16, 64, 256)
Mn = themap(Symplectic(order = 4, n_steps = n))
println("n_steps=", lpad(n,3), ": kx = ", -Mn[2,1], " ky = ", -Mn[4,3])
end
```
On `main` (BeamTracking 0.8.0, Beamlines 0.10.1):
| | `kx` | `ky` |
|---|---|---|
| analytic | -1.4908364 | 0.6618504 |
| `SciBmadStandard` (default) | -1.4974041 | 0.6821204 |
| `Symplectic(order=4, n_steps=1)` | -1.4974041 | 0.6821204 |
| `n_steps=2` | -1.4917217 | 0.6627882 |
| `n_steps=8` | -1.4908406 | 0.6618538 |
| `n_steps=64` | -1.4908364 | 0.6618504 |
The default matches `n_steps = 1` exactly, and is 3.1% high in `ky`. It converges to analytic by `n_steps ≈ 8-16`.
## Effect on a real lattice
The CESR lattice from bmad-sim/SciBmad.jl#87 has twelve of these magnets (`dqx1b/1d/2b/2d/3b/3d/4b/4d/5b/5d/6b/6d`), all identical. One-turn matrix on `main`, changing nothing but the step count on those twelve elements:
```
default (n_steps = 1): trace/2 x = -0.9337096 trace/2 y = +0.9808064 fractional tunes = [0.031233, 0.441724]
n_steps = 16: trace/2 x = -0.9386931 trace/2 y = -0.6515009 fractional tunes = [0.362930, 0.443981]
n_steps = 64: trace/2 x = -0.9386880 trace/2 y = -0.6515117 fractional tunes = [0.362932, 0.443979]
```
The ring is stable either way, so nothing errors — but the default reports a **vertical fractional tune of 0.031 instead of 0.363**. The horizontal tune is roughly right, which makes the result look credible.
Everything else in that lattice is fine: the 122 quadrupoles get exact thick-lens maps (e.g. `qx4a`, `Kn1=1.53946`, `L=0.4` tracks -0.641375 vs analytic -0.641373), and the three wigglers match the analytic planar-wiggler focusing `k_y = (B₀/Bρ)²/2` to 0.16% and are converged in step count. Only the combined-function bends are affected, because they are the only elements in the lattice with no closed-form map.
## On released v0.7.0 the same defect makes the lattice unstable
v0.7.0 — which SciBmad currently resolves to — uses a worse split for `bend_bquadrupole`, and one step gives:
```
0.7.0, n_steps = 1: kx = -0.6957835 ky = 1.1296044 (+71% in ky)
analytic: kx = -1.4908364 ky = 0.6618504
```
≈ +0.47 m⁻¹ of spurious vertical focusing per magnet, ×12, which drives the one-turn matrix to `trace/2 y = -19.73` and makes `twiss` fail with `At least one orbital eigenmode is linearly unstable!`. Two steps is already enough to restore stability. This is the second problem behind SciBmad.jl#87 (the first being the ForwardDiff partial-count crash, bmad-sim/BeamTracking.jl#326).
The converged tunes agree between patched-0.7.0 and `main` to five decimals (0.362937/0.443981 vs 0.362932/0.443979), so both versions integrate the same physics correctly once given enough steps.
## Why this is easy to miss
A symplectic integrator with too coarse a step yields a map that is wrong but *exactly* symplectic — `det = 1` to 12 digits in every case above. Symplecticity checks therefore pass, and on `main` the eigenvalues stay on the unit circle too, so the only symptom is a wrong number.
## Suggestion
Choose the fallback step count from the element's integrated strength rather than fixing it at 1 — e.g. enough steps to keep `√|k₁|·ds` (and the analogous quantity for higher multipoles) below ~0.1, with 1 step as the floor. Failing that, a warning when a single-step fallback is applied to an element whose integrated strength exceeds some threshold would at least make it visible.
Worth noting for translated lattices specifically: Bmad's `bmad_standard` gives combined-function magnets an analytic map, so a Bmad lattice whose element data translates perfectly can still track very differently in SciBmad, with no indication that anything is off.
## Environment
Julia 1.11.7, aarch64-apple-darwin. BeamTracking 0.8.0 (`main`) / 0.7.0, Beamlines 0.10.1 / 0.9.4, GTPSA 1.5.6.
Contributor guide
No contributing guide indexed for this repository
Research direction
Start at SciBmadStandard's fallback to Symplectic() and reproduce the combined-function SBend case using the Julia snippet, comparing n_steps values with the analytic kx and ky results. Done means the fallback no longer silently uses an inadequate single step for strong elements, the reproduced tunes converge to the analytic results, and a regression check covers the affected lattice behavior.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- julia
- Domain
- backend
- Issue type
- Bug
- Difficulty
- 4/5
- Estimated time
- 3-5 days
- Activity status
- Active
- Clarity
- Mostly clear
- Newbie friendliness
- 55/100