544 closed orbit finder for opalx - #553
Conversation
Expose validated map-only OPTION settings for six perturbation amplitudes, Richardson levels 0..4, and typed ray-integrator dispatch (Boris/LF2). Extrapolate centered segment maps before ordered composition and record effective settings and refinement diagnostics. Advance host reference and shadow rays directly in metres, integrate signed path from speed, retain compensated position/time/path state through sampling and coordinate differences, and localize start/stop path targets by tracked bisection. Keep the Boris kick, boundary tolerances, production particle pusher and GPU kernels unchanged. Preserve the full reproducible DBA DT/amplitude/level studies, stdout, input hashes, source snapshots, analytic all-36-entry comparisons, MPI repeats, plots and audit scripts. Across 385 configurations per baseline, small-amplitude robustness improves substantially but peak scaled max-entry accuracy changes from 7.48e-9 to 8.37e-9; this is not an accuracy-floor improvement. Defaults and numerical tolerances are unchanged. Validation: 87/87 CTests; 20/20 focused tests on one and two MPI ranks; 12/12 Python study tests; four analytic optics regressions; 13 identical MPI map/diagnostic repeats; Doxygen and rendered manual formulas checked. Unrelated user input changes are excluded. Manual updates are committed separately in opalx-manual.
Extend the existing map integrator selector without changing the production particle pusher. Integrate the relativistic Lorentz equations and path length at each RK stage, preserve compensated position/time/path updates, and resolve field-support transitions independently for every ray. Document the RK tableau, units, stage cost, fixed-DT behavior and lack of a symplecticity guarantee. Cover analytic relativistic helices and convergence order, time-dependent electric fields, reverse stepping, material hits, drift maps, support overlap and option isolation. Validation: OPALX and TestOrbitThreader build; TestBendFieldModel and TestOrbitThreader CTest targets pass with OMP1. All 25 OrbitThreader cases also pass in a dedicated two-rank MPI run. These integrators are prerequisites for the recorded map-3 DBA comparison.
Publish the map-3 DBA input with HGAP=0.1 m, FINT=0.1, separated fringe supports and the existing unmatched quadrupole seed. Include analytic field profiles, all three single-rank integrator cases with stdout/timing/maps, source and artifact provenance, and the full 6x6 variational reference. Separate distributed-field numerical errors from the CERN thin-edge approximation and document the current full-gap/half-gap convention mismatch without modifying fringe physics. Retain the original analysis scripts for historical hash verification; default replay uses the published settings without requiring private map-2 campaign paths. Add the shared field-history parser and configurable map-options helper needed by the study. Keep intermediate experiments, HDF5 output, executables and unrelated sandbox changes out of the commit. Link the authoritative physics manual. Validation: all 25 Python map-3 tests pass both in the worktree and an isolated export of the staged files, plus 4 SDDS parser tests. Recomputed full-map differences reproduce the recorded comparison; no numerical tolerances changed. Automated OPALX regression integration remains follow-up work.
Keep nominal lattice geometry, fringe-field support and travelled reference path length distinct. Fringe fields may contribute inside a drift without extending its nominal length or changing map ownership. - Add nominal-body containment queries and partition transfer maps by body ownership, while retaining summed-field evaluation for all shadow rays. - Record field contributors separately from map owners, including fringe fields in drifts and unowned intervals. - Populate IndexMap from support-resolved reference steps and coalesce adjacent identical support sets without merging separate visits. - Calculate RING design circumference recursively from nominal occurrence arc lengths, including drifts and repeated nested cells. - Locate the same-direction return to the starting transverse plane using a bounded search, rather than stopping at the design circumference. - Use the measured return length for IndexMap periodicity and reuse the entrance axes for the return map. - Report design circumference, return length, position displacement and relative momentum mismatch. - Resolve short nominal bodies in ray stepping and reject premature launch-plane crossings during map construction. - Update Doxygen descriptions of ownership, field contributions and same-section return-map semantics. Add unit coverage for circumference independence from fringe extent, fringe contributions inside nominal drifts, cutoff convergence, circular reference returns, missing-return failures and IndexMap coalescing. This does not implement closed-orbit finding or change the native FINT model. The independently observed finite-fringe DBA map discrepancy remains unresolved. Refs #546
Remove all 3345 sandbox paths from the repository index. These files were additions on this branch and are not present on master. Keep source, unit tests and other tracked files unchanged. Ignore /sandbox/ to prevent local study inputs, generated output and exploratory work from being staged accidentally again. Use index-only removal so local files and uncommitted edits remain intact. Verify the complete local sandbox path inventory and all 8927 file content hashes are unchanged. This is a normal cleanup commit, not a history rewrite: earlier commits still retain their sandbox snapshots.
Extend RFCAVITY TYPE=SINGLEGAP with a separate radial voltage-profile model, leaving continuous-field cavity tracking unchanged. Add strict profile loading, cubic Hermite interpolation using supplied derivatives, and shared host/device kick formulas for transit-time energy gain and legacy focusing. Document SI geometry, phase conventions, and the retained legacy extra 2*pi focusing factor. Add a deliberately scoped one-proton, one-container, one-rank event path in ParallelTracker. Resolve directed offset-gap crossings with 40 Boris-substep bisections, apply full impulses in chronological order, and track the remainder. Use Cartesian element poses and spatial sector support rather than the nominal closed-orbit index. Transform the particle's moving coordinates to/from the lab frame and use the same integration routine for particle and reference. Reject unsupported species/spin/restart/field-solver configurations. This is a host event path, not a production multi-particle or GPU event implementation; particle E/B diagnostic arrays are not populated. Keep CyclotronSector concrete without adding a Rep class, consistent with the planned removal of Rep classes. Provide isolated old-OPAL and OPALX first-turn generators/comparisons, original accelerating input and voltage profiles, and compact validation summaries. Original input uses obsolete inline trim syntax; adaptation follows the named mirrored coil in cyclotron1.in. Historical trim equivalence remains to be checked before relying on high-radius TC comparisons. Validation (OpenMP CPU, one MPI rank; build -j 10): - 72 -> 73.7423153 MeV with five expected RF kicks at DT/4 (10.283 ps). - Maximum old-LF2 position difference decreases 0.758 -> 0.507 -> 0.186 um. - Finest-step final energy difference -1.93 eV; particle/reference <12.1 pm. - TestCyclotronRF (six cases), TestRFCavity, TestCyclotronSector and TestOrbitThreader pass; full coasting study retains 8.53 nm LF2 agreement. - Kernel Doxygen validation and separately updated manual chapter renders pass. Energy-target stopping and full acceleration to 590 MeV are the next milestone. GPU execution, multi-particle event tracking and restart are not validated. Refs #545
… 590 MeV Extend TRACK with EKINSTOP in GeV, transferred internally in eV. Stop at the first complete reference RF kick reaching the requested kinetic energy; retain the physical overshoot and omit the remaining magnetic drift. Precompute the reference terminal substep before advancing the bunch clock, then publish it through the normal reference/frame update. Advance the physical particle independently rather than substituting the reference state. Require the non-restarted single-proton SINGLEGAP RING path, one positive DT segment, a target above launch, and no explicit RUN TURNS. Keep MAXSTEPS a hard safety budget and reject exhaustion or loss before target. Continue reporting directed return-plane turns; suppress nominal-circumference progress for gaps. The default disabled target preserves existing tracking behavior. Extend reproducible one-rank old-OPAL/OPALX runners for full acceleration, TC-off/on and matched fixed-step endpoints. Explicitly set old OPAL trim-coil PHIMIN=0, PHIMAX=360: the installed 2022.1 constructor has no PHIMAX default despite its help text. Preserve original inputs and do not modify old source. 72 -> 590 MeV results, 949 OPALX gap events in each endpoint comparison: - No TC, 41.1319513 ps: old 590.325895173, new 590.326824538 MeV; +929.365 eV, 55.785 um final position difference. - No TC, half DT: old 590.296910517, new 590.297352913 MeV; +442.396 eV, 26.682 um final position difference. - Active TC, full DT: old 590.385013924, new 590.386116234 MeV; +1102.310 eV, 63.183 um final position difference. - Particle/reference separation remains below 3.2 nm. Compare coordinates at matched step endpoints, not the earlier EKINSTOP gap time. Report empirical inter-code guards (1.5 keV, 100 um) separately from absolute accuracy: old LF2 itself shifts 29 keV on halving DT. No asymptotic order or active-TC convergence claim is made. Retain compact summaries and sampled errors/energies; exclude large raw trajectories and generated plots. Validation: build -j10; TestCyclotronSector, TestCyclotronRF, TestRFCavity and TestOrbitThreader all pass; seven energy-stop rejection cases; positive short and full energy-target runs; first accelerating turn unchanged at 0.758 um; coasting DT refinement rerun with old-LF2 agreement 8.5256 nm; report guards pass; RF-kernel Doxygen warnings-as-errors and git diff --check pass. Expand Doxygen for cavity modes, units, stopping semantics and pending state. User and physics chapters were also updated in the separate local manual repo, validated and rendered (not published by this source commit). Remaining scope: one unpolarized median-plane proton, one rank/container, no restart or production GPU event path, E/B diagnostic arrays not populated. No new Rep subclass, upstream dependency changes or numerical tolerance changes to existing tests. Follow-up: active-TC timestep refinement/generalization. Refs #545
Theory and numerical method
---------------------------
- Introduce a fixed-energy, four-dimensional closed-orbit finder solving
T(u) - u = 0, with u = (x, px, y, py).
- Express positions in metres and mechanical momenta in units of mc.
- Define one turn by a directed return to a fixed transverse section,
allowing each perturbed trajectory its own return time.
- Construct the one-turn Jacobian with central finite differences.
- Use a scaled, damped Newton iteration and pivoted QR solve, requiring
both residual and correction convergence plus a fresh verification turn.
- Report the closed orbit, one-turn matrix, eigenvalues and fractional
mode tunes to stdout and result files.
LAPACK eigenanalysis
--------------------
- Fetch LAPACK/LAPACKE through CMake FetchContent.
- Use LAPACKE_dgeev for complex eigenpairs of the real, nonsymmetric
transverse matrix, with backend error and eigenvector-residual checks.
- Use singular-value analysis to diagnose an ill-conditioned eigenbasis.
- Distinguish stable, unstable, non-unit-circle and marginal spectra.
- Report conjugate fractional phases without inventing integer tune
branches or assigning coupled modes to transverse planes.
- Retain the current Fortran-backed LAPACK implementation; record a
possible future C-only backend as follow-up work.
COF command and ring-aware aperture checking
-------------------------------------------
- Add the scoped COF / RUN / ENDCOF input interface for static magnetic
RING configurations, without particle emission or collective tracking.
- Preserve native element fields and shared ray integrators.
- Prevent false aperture losses from opposite ring-arm longitudinal
slabs by selecting the nearest finite conventional design centreline,
independently of aperture dimensions.
- Retain genuine native aperture-loss checks and cyclotron-sector domain
checks; document the local-path assumption and intersecting-pipe limit.
Stable, non-achromatic DBA-ring benchmark
-----------------------------------------
- Add sandbox/dba-ring-cof, preserving historical achromatic examples.
- Generate a six-cell RING from explicitly placed elements using
X/Y/Z/THETA/PHI/PSI.
- Each cell contains:
SBEND: radius 2 m, angle 30 degrees, HGAP=0, FINT=0
DRIFT: length 1 m
QUADRUPOLE: length 0.2 m, K1=+2.84 m^-2
DRIFT: length 1 m
SBEND: radius 2 m, angle 30 degrees, HGAP=0, FINT=0
- Use electrons with PC=0.2505104781131461 GeV/c.
- Preserve the 25.76637061436 m design circumference.
- Select stable optics analytically without imposing achromaticity or
introducing a second quadrupole family.
- Provide generated COF inputs, analytic scans, comparison scripts,
saved matrices, convergence studies and run provenance.
Analytic references and coordinate conventions
----------------------------------------------
- Reuse the exact hard-edge sector-bend, drift and thick-quadrupole
transfer matrices in sandbox/map-2/check_maps.py.
- Form the analytic one-turn reference as M_ring = M_cell^6.
- Compare like coordinates using M_u = T M_slope T^-1, where
T = diag(1, p0/(mc), 1, p0/(mc)), evaluated with the actual BEAM momentum.
- Keep this exact hard-edge reference distinct from the distributed-field
numerical reference and thin-edge approximation in sandbox/map-3/dba.
DBA validation results
-----------------------
- Run physics studies with one MPI rank and one OpenMP thread.
- With DOP853, DT=5e-12 s and central-difference position/slope amplitudes
of 3e-5, obtain:
Analytic Tracked
Horizontal 0.229892430867 0.229892431493
Vertical 0.310463602980 0.310463603175
- Maximum scaled transverse matrix-entry error: 6.32e-9.
- Maximum eigenvalue-modulus error: 5.13e-10.
- Both modes are STABLE under the unchanged 1e-8 unit-circle tolerance.
- Confirm agreement at DT=1e-11 s and convergence from a displaced seed.
- Retain finite-difference scans showing increased numerical noise at
smaller perturbations; do not relax tolerances to force stability.
Tests, documentation and remaining scope
----------------------------------------
- Pass seven selected CTest suites and three analytic benchmark tests.
- Cover opposite-arm false losses, genuine local aperture losses,
parser errors, repeated COF blocks and dedicated two-rank consistency
and failure propagation.
- Update source documentation and user/physics manual descriptions.
- Keep dispersion, chromaticity, integer tune reconstruction, spectral
DBA cross-checks and finite-fringe ring validation as follow-up work.
Replace the GoogleTest-scoped Kokkos lambda with a namespace-scope functor accepted by NVCC. Limit project warning flags to C and C++ so fetched Reference LAPACK Fortran sources do not inherit them, and suppress the external LAPACKE C-linkage warning locally. Refs #544
Expose the cost of reference tracking, segmentation, map construction and output, then reduce repeated membership searches and private-ray exit localization work. Add phase timers and counters for field evaluations, element lookups, recursive subdivisions, accepted steps, membership reuse and exit trials. Report executable identity, build configuration, compiler flags, assertions, MPI ranks and OpenMP/Kokkos concurrency to make timing comparisons auditable. Correct build-type lookup when recording compiler optimization flags. Reuse element membership only at structurally identical recursive states. Preserve lazy lookups, field evaluation order and the accepted integration tree. Regression tests compare cached and uncached trajectories bit for bit across integrators, overlapping supports and forward/reverse tracking. Replace fixed exit-plane bisection with safeguarded secant proposals for RK4 and DOP853. Reintegrate every trial from the original bracket start, force bisection after poor contraction, and retain the 1e-12*abs(DT) time tolerance. Reject invalid or unconverged brackets. Keep the original bisection sequence for Boris, including its LF2 alias. The stable DBA study reproduces material Boris sensitivity to changed exit trial durations. Neither search is established as the accuracy winner. The final Boris fallback reproduces the original map at printed precision and retains all 14,400 exit iterations in the matched DBA check. Add the stable hard-edge DBA reference, reproducible study scripts, inputs, compact results and executable/input provenance. Include timestep and finite-difference checks for Boris, RK4 and DOP853, a displaced-start COF check, and comparisons with analytic optics and direct COF Jacobians. Exclude executables, full trajectories and unrelated sandbox outputs. Validation: - Release build and six relevant CTest suites pass. - Three independent analytic-reference tests pass. - All 15 COF grid runs verify closure; displaced-start recovery also passes. - At DT=5e-12 s and h=3e-5, maximum tune errors are 1.22e-6 for Boris, 4.56e-10 for RK4 and 6.26e-10 for DOP853. - Stable-DBA exit time falls from 7.17 to 0.433 s for RK4 and from 20.28 to 1.39 s for DOP853 in individual sequential runs. - Old/new reference trajectories agree exactly in the matched studies. No public API, integration method, numerical tolerance or parallel execution order is changed. Secant trial times can change final roundoff. Timing samples are not statistical speedup estimates, and the hard-edge DBA results do not establish accuracy for the ISIS lattice.
Replace the temporary initializer-list loop over FDSTEP and SCALES with a named std::array. Preserve option defaults, array-length validation and value assignment without suppressing compiler warnings. CDash build 4048748 reported two -Wdangling-reference warnings and their accompanying notes at this loop. Validation: - GCC 15 reproduces both warnings with the original loop. - The replacement passes with -Werror=dangling-reference. - The local Release build passes. - The single-rank, four-thread DBA COF result is exactly unchanged, including coordinates, residuals, matrix, eigenvalues and tunes. The remote CDash build has not yet been rerun.
ute help
The space-charge mesh and particle positions are transformed into the beam-aligned coordinate system before the field solve. However, particle momenta previously remained in the reference frame. This mixed coordinate systems when computing bin membership, mean bin momentum, relativistic gamma factors, and the Lorentz transformation of the self-fields. The inconsistency is particularly relevant for beams following curved reference trajectories. Rotate particle momenta into the charge-mesh frame before redistribution, binning, and the self-field solve. Rotate them back into the reference frame after the solve, together with the existing position and field transformations. Apply the momentum transformation only when binning is active, leaving the legacy non-binned path unchanged. The rotations execute on the device and introduce no additional host transfers or MPI collectives. They preserve momentum magnitude and gamma up to floating-point roundoff. Document that binned self-field calculations require positions and momenta to use the same coordinate frame. Add a regression test covering straight and rotated beam frames, zero-momentum particles, momentum-norm preservation, and the complete rotation round trip. A one-turn 72 MeV cyclotron benchmark using 262144 particles, a 32^3 mesh, one momentum bin, and the STANDARD Green function gives final OPALX-to-OPAL differences of approximately 0.07% or less for all RMS sizes, 0.17% or less for transverse emittances, 0.67% for longitudinal emittance, and 0.41% for RMS energy spread. Detail can be found #552
Add the general RUN attribute SCFIELDUPDATE to control when
space-charge fields are evaluated within the drift-kick-drift
integration step.
Support two update modes:
- MIDPOINT evaluates the self-field after the first half drift, using
particle positions at R_(n+1/2). This remains the default and
preserves existing OPALX behavior.
- PRESTEP evaluates the self-field before the first half drift, using
particle positions at R_n. The gathered particle fields are carried
through the half drift and combined with external fields evaluated at
R_(n+1/2), reproducing the historical OPAL space-charge ordering.
The Boris momentum update, position drifts, particle deposition, field
solver, Green function, and external-field evaluation remain unchanged.
SCFIELDUPDATE is a general parallel-bunch tracking option and is not
specific to rings or cyclotrons. It has no physical effect when the
selected field solver disables space charge.
Example:
RUN, METHOD="PARALLEL", FIELDSOLVER=FS0,
SCFIELDUPDATE="PRESTEP";
Validate PRESTEP with a one-turn 72 MeV cyclotron benchmark using a
2 mA beam, 262144 particles, a 32^3 mesh, one momentum bin, the
STANDARD Green function, one MPI rank, and one OpenMP thread.
Relative to old OPAL, changing from MIDPOINT to PRESTEP changes the
final differences as follows:
- radial RMS size: +0.0024% to -0.0041%
- radial emittance: -0.3118% to -0.2984%
- vertical emittance: +0.1167% to -0.0153%
- longitudinal emittance: -1.3186% to -0.8755%
- RMS energy spread: -1.4927% to -1.0705%
The comparison shows that historical field timing accounts for a
significant part of the vertical and longitudinal differences, while
remaining discrepancies require separate investigation. Again the results
are documented here #552
|
I checked a bit the branch In the more complete branch However, here, it seems that the tracking stops at the first integration endpoint after the return-plane crossing, while |
|
In case that this was not intentional, the space charge examples that were attached on Element run for 720 time steps and this seems to be a bit too short to complete the full turn. |
|
which one? |
Integrate the space-charge solver refactor while retaining ring controls and PRESTEP/MIDPOINT timing. Adapt the COF reference container and preserve SINGLEGAP solver validation in RUN.
Locate the final directed return-plane crossing using TRACK’s Boris reference step and advance the bunch with a shared shortened timestep. Keep MPI ranks, particle times and final output consistent. Add unit and MPI regression tests and ISIS timestep studies. Validate COF-to-TRACK initialisation: the transverse closure residual drops from 512 µm to 3.21 µm at 100 ps and to 0.185 µm at 25 ps. Limit event localisation to static magnetic, single-container runs without space charge.
Both the midpoint and the prestep examples. |
|
With the latest commits, the overshoot from the last time step seems to be fixed. Also for my previous comment, adding |
Purpose
General fixed-energy RING closed-orbit finding, with a one-turn matrix and complex-eigenvalue/tune diagnostics. Native magnetic rings and cyclotron sectors share the same solver.
Related: #544 and #545. This branch also contains preceding cyclotron, spectral-tune and transfer-map work relative to master; the focused COF commit is b8de53d.
Regression tests
The Ring and DBA test will be the regression tests after this is merged.
Theory
Use section coordinates$u=(x,p_x,y,p_y)^T$ , with positions in metres and mechanical momenta $p_i=P_i/(mc)$ . At fixed total momentum $p_0$ :
with scaled, damped Newton steps:
Column-pivoted QR solves the Newton system. Both residual and undamped correction must satisfy position/momentum tolerances (default$10^{-10}$ in the respective units), followed by a fresh verification turn.
LAPACK/LAPACKE: use dgeev for the full real nonsymmetric matrix's complex eigenpairs and dgesvd to diagnose eigenbasis conditioning. Check backend status, eigenvector residuals, conjugate pairing and distance from the unit circle. Stable modes yield fractional phases and their complements; integer branches and automatic transverse-plane labels are not inferred for general coupled maps.
Reference LAPACK/LAPACKE 3.12.1 is checksum-pinned and fetched with CMake. A Fortran compiler is currently required; no separate LAPACK installation is needed. Keep this checked backend until IPPL supplies Kokkos Kernels together with the required numerical backend.
This is a 4D fixed-momentum solve, not 6D synchronous closure. Dispersion is not assumed zero: momentum dependence is simply outside the current four-variable solve.
Input sketches
The snippets deliberately omit repeated placements; they illustrate the input structure rather than constitute complete runnable decks.
PSI Ring: 72 MeV coasting proton
No RF, particle emission or collective solve. The quoted unused source satisfies BEAM syntax without loading a distribution. The default return section is the RING frame; an explicit SECTION={X,Y,Z,THETA,PHI,PSI} is supported (metres/radians).
A separate compact native-ring regression uses four 90-degree SBENDs: radius 1 m, L=PI/2, ANGLE=PI/2, K1=-0.3, successive ELEMEDGE positions, and R: RING=(B0,B1,B2,B3). Proton and electron cases are tested.
Stable, non-achromatic DBA ring
Six cells, each B1–D1–Q–D2–B2:
Momentum perturbations correspond to$p_0,3\times10^{-5}$ , i.e. a slope-equivalent amplitude, not a mechanical-momentum amplitude of 3e-5.
Analytic reference
The original achromat setting K1=-6.371966681365967 is transversely unstable. The new setting preserves geometry but deliberately allows dispersion.
Use exact hard-edge drift, sector-bend and thick-quadrupole matrices:
For focusing$k>0$ ,
Use its hyperbolic continuation for$k<0$ and drift limit for $k=0$ . In this electron convention $Q_x=Q(-K_1,L_Q)$ and $Q_y=Q(K_1,L_Q)$ .
Compose in traversal order and calculate$M_{\rm ring}=M_{\rm cell}^{6}$ . For the on-axis orbit:
The conversion uses the actual BEAM normalised momentum.
Reference implementations: local sandbox/map-2/check_maps.py (drift, sector_bend, quadrupole, dba), and sandbox/dba-ring-cof/analytic_scan.py and run_benchmark.py. These sandbox files are not in the currently pushed tree; publishing reproducible fixtures is an open item. The equations and numerical comparison are therefore included directly here. The distributed-field reference and thin-edge approximation under local sandbox/map-3/dba are different models, not exact references for this hard-edge ring.
Results: DBA map and tunes
One MPI rank, one OpenMP thread, DOP853, DT=5e-12 s, central-difference position/slope amplitude 3e-5.
Both matrices are in (x,x',y,y'). All cross-plane entries are zero for this on-axis median-plane case.
M12/M34 have units m; M21/M43 have units 1/m. With positions scaled by 1 m, maximum dimensionless matrix-entry error is 6.32e-9.
Results: PSI Ring command example
At the coarser 72 MeV RK4 settings above:
The mechanical-coordinate map is approximately:
Orbit closure succeeds, but the full spectrum is NON_UNIT_CIRCLE at these settings. This is not a validated stable radial tune. Cyclotron timestep/FD sensitivity remains distinct from the successful hard-edge DBA validation.
Correctness and validation
Open points
Dispersion, chromaticity, finite-fringe DBA and general 6D validation are not claimed here.