Generalize 6D canonical integrator to arbitrary curvilinear coordinates (VMEC)

CMakeLists.txt

@@ -12,7 +12,12 @@ enable_language(C Fortran)
 # Optional: enable CGAL-based STL wall intersection support.
 option(SIMPLE_ENABLE_CGAL "Enable CGAL STL wall intersection support" OFF)
 
-# Get version from git describe (format: v1.5.2-42-g62fed6e[-dirty])
+# Version from git describe (format: v1.5.2-42-g62fed6e[-dirty]). Generated at
+# configure time so the file exists for the first build, AND regenerated at BUILD
+# time by an always-run target (cmake/GenerateVersion.cmake) so the baked-in
+# string follows the current HEAD/dirty state instead of going stale after a
+# commit without a reconfigure. The script rewrites version.f90 only when the
+# describe output actually changes, so an unchanged HEAD triggers no recompile.
 execute_process(
     COMMAND git describe --tags --dirty --always
     WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}
@@ -25,13 +30,23 @@ if(NOT SIMPLE_VERSION)
 endif()
 message(STATUS "SIMPLE version: ${SIMPLE_VERSION}")
 
-# Generate version.f90 from template
 configure_file(
     ${CMAKE_SOURCE_DIR}/src/version.f90.in
     ${CMAKE_BINARY_DIR}/version.f90
     @ONLY
 )
 
+add_custom_target(simple_version_gen ALL
+    COMMAND ${CMAKE_COMMAND}
+        -DSRC=${CMAKE_SOURCE_DIR}/src/version.f90.in
+        -DDST=${CMAKE_BINARY_DIR}/version.f90
+        -DGIT_DIR=${CMAKE_SOURCE_DIR}
+        -P ${CMAKE_SOURCE_DIR}/cmake/GenerateVersion.cmake
+    BYPRODUCTS ${CMAKE_BINARY_DIR}/version.f90
+    COMMENT "Refreshing SIMPLE version from git describe"
+    VERBATIM
+)
+
 add_compile_options(-g)
 
 # Compiler-specific flags
@@ -164,6 +179,15 @@ else()
     find_package(LAPACK REQUIRED)
 endif()
 
+# The 6D full-orbit / CPP integrators need the metric tensor + Christoffel
+# symbols from libneo (issue #322). Until that merges to libneo main, pin the
+# fetched libneo to the feature branch by default. Override with -DLIBNEO_REF=...
+# (e.g. main once merged) or -DLIBNEO_SOURCE_DIR=... for a local checkout.
+if(NOT DEFINED LIBNEO_REF AND NOT DEFINED LIBNEO_SOURCE_DIR)
+    set(LIBNEO_REF "feature/metric-christoffel" CACHE STRING
+        "libneo git ref to fetch (branch, tag, or SHA)")
+endif()
+
 find_or_fetch(libneo)
 
 # Consume the NetCDF stack libneo resolved (see comment near the top). In the
@@ -439,6 +463,10 @@ add_executable(diag_traj.x
     app/simple_diag_traj.f90
 )
 
+add_executable(diag_cp_gc.x
+    app/simple_diag_cp_gc.f90
+)
+
 # Apply SIMPLE-specific compile options to executable
 target_compile_options(simple.x PRIVATE ${SIMPLE_COMPILE_OPTIONS})
 
@@ -459,6 +487,7 @@ target_compile_options(diag_meiss.x PRIVATE ${SIMPLE_COMPILE_OPTIONS})
 target_compile_options(diag_albert.x PRIVATE ${SIMPLE_COMPILE_OPTIONS})
 target_compile_options(diag_orbit.x PRIVATE ${SIMPLE_COMPILE_OPTIONS})
 target_compile_options(diag_traj.x PRIVATE ${SIMPLE_COMPILE_OPTIONS})
+target_compile_options(diag_cp_gc.x PRIVATE ${SIMPLE_COMPILE_OPTIONS})
 
 # Apply trampoline error flags only to SIMPLE executable (not subprojects like fortplot)
 if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
@@ -481,12 +510,28 @@ target_link_libraries(diag_albert.x PRIVATE simple)
 target_link_libraries(diag_newton.x PRIVATE simple)
 target_link_libraries(diag_orbit.x PRIVATE simple)
 target_link_libraries(diag_traj.x PRIVATE simple)
+target_link_libraries(diag_cp_gc.x PRIVATE simple)
 set_target_properties(diag_traj.x PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
+set_target_properties(diag_cp_gc.x PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
 set_target_properties(diag_meiss.x PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
 set_target_properties(diag_albert.x PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
 set_target_properties(diag_newton.x PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
 set_target_properties(diag_orbit.x PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
 
+# Boozer chartmap export tool
+add_executable(export_boozer_chartmap.x
+    app/export_boozer_chartmap.f90
+)
+target_compile_options(export_boozer_chartmap.x PRIVATE ${SIMPLE_COMPILE_OPTIONS})
+if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
+    target_compile_options(export_boozer_chartmap.x PRIVATE
+        $<$<COMPILE_LANGUAGE:Fortran>:-Wtrampolines>
+        $<$<COMPILE_LANGUAGE:Fortran>:-Werror=trampolines>
+    )
+endif()
+target_link_libraries(export_boozer_chartmap.x PRIVATE simple)
+set_target_properties(export_boozer_chartmap.x PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
+
 # Ensure fortplot is built before diagnostics (when available)
 if(NOT CMAKE_Fortran_COMPILER_ID STREQUAL "NVHPC")
     add_dependencies(diag_meiss.x fortplot)

DOC/coordinates-and-fields.md

@@ -586,6 +586,224 @@ f%dAth = [Ath_norm, 0, 0] ! Constant derivative
 - Uses libneo functions: `vmec_to_can`, `can_to_vmec`
 - Simpler than Meiss/Albert but less optimized
 
+### 6.6 The Curvilinear 6D Canonical Integrator
+
+**Files**: `src/orbit_cpp_canonical.f90`, `src/orbit_cpp_vmec_metric.f90`,
+`src/orbit_cpp_chartmap_metric.f90`, `src/field/field_can_test.f90`
+(`eval_field_correct_test`)
+
+The guiding-center integrators reduce the perpendicular motion to the magnetic
+moment. The 6D canonical integrator in `orbit_cpp_canonical` keeps the full
+phase space `(q, p)` and resolves (or, for the Pauli models, represents) that
+motion directly. It is the SIMPLE port of the Egger-Feiel thesis
+discrete-variational integrators, generalized to arbitrary curvilinear
+coordinates with a full (non-diagonal) metric.
+
+The Hamiltonian is `H = (1/2m)(p_i - qc A_i) g^ij (p_j - qc A_j) [+ mu|B|]`, so
+`q_dot^k = (1/m) g^kj (p_j - qc A_j)` and
+`p_dot_k = qc A_{j,k} v^j + (m/2) g_{ij,k} v^i v^j [- mu |B|_{,k}]`. Every term
+carries the full metric `g_ij`, its inverse `g^ij`, and the direction
+derivatives `g_{ij,k}`. The integrator reads them from a `block_t`: metric,
+metric derivatives, covariant `A_i` with gradient, `|B|` with gradient, and the
+covariant unit field `h_i`.
+
+Three coordinate blocks fill that structure.
+
+`COORD_TOK` is the analytic tokamak, inline and GPU-portable. The metric is
+diagonal, `g = diag(1, r^2, (R0 + r cos theta)^2)`,
+`sqrt(g) = r (R0 + r cos theta)`. The covariant vector potential is `A_r = 0`,
+`A_theta = B0 (r^2/2 - r^3 cos(theta)/(3 R0))`,
+`A_phi = -B0 iota0 (r^2/2 - r^4/(4 a^2))`, with `B0 = iota0 = R0 = 1`, `a = 0.5`.
+The 6D path needs the exact field from the curl of `A`:
+`B^k = eps^ijk A_{j,i} / sqrt(g)`, `|B| = sqrt(g_ij B^i B^j)`, so with `A_r = 0`,
+`|B|^2 = A_{phi,r}^2 / (R0 + r cos theta)^2 + A_{theta,r}^2 / r^2`
+(`eval_field_correct_test`). The guiding-center path keeps the linearized
+`eval_field_test` (`|B| = B0(1 - r/R0 cos theta)`); at the reference start
+`(r,theta) = (0.1, 1.5)` the exact `|B| = 0.99749` against the linearized
+`0.99293`. The diagonal metric is the special case of the general arithmetic
+(off-diagonals zero), so `COORD_TOK` reproduces the python oracle bit-for-bit
+while the same residual runs on a stellarator metric.
+
+`COORD_VMEC` runs on real VMEC equilibria in native flux coordinates
+`(s, vartheta, varphi)`, wired through `orbit_cpp_vmec_metric`. It is the
+production 6D loss chart (see below). The full metric `g_ij`, `g^ij` and
+Christoffel symbols `Gamma^l_jk` come from libneo's `coordinate_system_t`
+(issue #322, branch `feature/metric-christoffel`); the metric derivatives follow
+from metric compatibility, `g_{ij,k} = g_il Gamma^l_jk + g_jl Gamma^l_ik`. The
+covariant `A_i` and `|B|` come from SIMPLE's native VMEC field
+(`vmec_field_evaluate`), with `dA` and `d|B|` by central difference. The metric
+is consistent with the field: the covariant unit field obeys
+`h_i g^ij h_j = |h|^2 = 1` to about `1e-2` (e.g. `0.998`, `1.008`, `0.946` at
+`s = 0.3, 0.5, 0.7`). The libneo metric inverts exactly, `|g g^-1 - I| ~ 6e-16`,
+so the residual is not a metric defect. It comes from the source mismatch: `h_i`
+is built from SIMPLE's native VMEC field `B_i`, while `g^ij` is libneo's
+independent metric, and the two splined representations of the same equilibrium
+differ at the percent level. (On the broken chartmap the same invariant was
+`O(nfp^2)`, hundreds, not `1.009`.)
+This block is host-side: libneo's metric is `class()`-dispatched and reads 3D
+splines, so it cannot run under `!$acc routine seq`.
+
+`COORD_CHARTMAP` runs the 6D state in `(rho, theta_B, phi_B)` with
+`rho = sqrt(s)`, wired through `orbit_cpp_chartmap_metric`: the chartmap metric
+and Christoffel from `reference_coordinates%ref_coords`, the field reparametrized
+from `s = rho^2` with the radial chain rule `dF/drho = 2 rho dF/ds`, and the
+covariant `A_i`, `h_i`, `|B|` from the active `field_can_mod%evaluate` pointer.
+This chart is NOT consistent and is not the production route. libneo splines the
+chartmap Cartesian `x/y/z` with a periodic fit over one field period, but for
+`nfp > 1` the Cartesian `x,y` are not field-period-periodic (they rotate by
+`2pi/nfp`), so the periodic spline destroys the secular toroidal rotation: the
+analytic spline `e_phi` loses its `~R` magnitude and the geometric metric gives
+`h_i g^ij h_j ~ 228..472` (`O(10^2)`) instead of 1. The defect is upstream in
+libneo's
+Cartesian-storage path and in the storage convention itself, so it cannot be
+repaired in the SIMPLE metric post-processor. A consistent chartmap route needs
+an R,Z (cylindrical) Boozer-chartmap representation in libneo: R and Z are
+field-period-periodic, and the reader's `has_spl_rz` path already adds the
+toroidal rotation analytically. Until then the production 6D loss path uses
+`COORD_VMEC`.
+
+The magnetic coupling carries a length normalization. The Hamiltonian uses
+`qc = charge/(c rho0)` with `rho0 = 1` by default, so `COORD_TOK` keeps the
+thesis `qc = charge/c` (`c = 1`). The production wire threads
+`rho0 = ro0_bar = ro0/sqrt(2)` so `qc = sqrt(2)/ro0`, which makes the canonical
+momentum `p_i = vpar h_i + A_i/ro0_bar` match the guiding-center `pphi` seed of
+`init_sympl`.
+
+Three models share one integer-dispatched residual: `MODEL_CP` (full charged
+particle), `MODEL_CPP_SYM` (Pauli symplectic midpoint, `H + mu|B|`),
+`MODEL_CPP_VAR` (Pauli variational midpoint, discrete Euler-Lagrange). `MODEL_CP`
+omits the `mu|B|` term: the full charged particle resolves the perpendicular
+gyromotion directly, so its kinetic energy already holds the perpendicular
+energy; the Pauli models drop the resolved gyromotion and reinstate it as the
+guiding-center `mu|B|`. The state is fixed-size 6,
+`z = (q1, q2, q3, p1, p2, p3)`; the position rows solve the canonical midpoint
+and the momentum rows carry `p`, so the Jacobian is square `6x6` and solved with
+the device LU `rk_solve` from `orbit_rk_core`. For `COORD_TOK` Newton uses the
+analytic Jacobian: `d2g` and `d2A` come from central differences of the block's
+own `dg`/`dA`, and the `O(mu)` `|B|` force gradient uses the block's analytic
+Hessian `d2Bmod`, the closed-form second derivative of the corrected `|B|`. The
+analytic-vs-finite-difference self-check passes for all three models. For
+`COORD_VMEC` the Jacobian is a central difference of the whole residual,
+consistent with the spline-based block; the FD step is taken relative to each
+variable's own scale (angles `O(1)`, momenta their own magnitude) so the
+physical-CGS momenta `~1e-8` keep an accurate column. A central-difference
+Jacobian is accurate to `~1e-7`, so the Newton step cannot reach the
+analytic-path floor `rtol = 1e-12`; the `COORD_VMEC` path converges on a
+step tolerance `rtol_fd = 1e-8` while `COORD_TOK`/`COORD_CHARTMAP` keep the
+tight `rtol`.
+
+The field `|B|` and its derivatives are the exact closed form of
+`|B| = sqrt(W)`, `W = A_phi,r^2/(R0 + r cos theta)^2 + A_theta,r^2/r^2`:
+`d_k|B| = dW_k/(2|B|)`, `d2_kl|B| = d2W_kl/(2|B|) - dW_k dW_l/(4|B|^3)`,
+finite-difference verified to `~1e-9` against `|B|`. An earlier port carried over
+two errata from the python listing. The metric theta-derivative now carries the
+factor `r`: `d g_33/d theta = -2 r (R0 + r cos theta) sin theta`. The field
+`d|B|/d theta` now keeps the `A_theta,rth` chain-rule term the listing dropped.
+With both corrected, the `CPP-sym` energy oscillation over a fixed time window
+converges as `dt^2`: `max|dE/E0| = 4.2e-5, 1.0e-5, 2.2e-6, 4.8e-7` for
+`dt = 80, 40, 20, 10`, each halving reducing the bound by about four. The buggy
+`d|B|` instead produced a flat `~1e-3` plateau that did not converge; the test
+asserts the `dt^2` behavior. The Fortran reproduces the regenerated python
+oracle for all three models to solver tolerance.
+
+`COORD_TOK` runs on the OpenACC device. `cpp_canon_step_tok` is the device entry:
+the whole chain (`residual_tok` -> `eval_block_tok` /
+`eval_field_correct_test` / `dLdq` / `raise` / `residual_blk`,
+`jacobian_analytic` -> `grad_jacobian_tok`, `rk_solve`) is `!$acc routine seq`
+with fixed-size state, integer model dispatch, the analytic Jacobian, and no
+`class()` or procedure pointer, so one particle runs per GPU thread. The host
+`cpp_canon_step` keeps the coordinate dispatcher; `COORD_VMEC` is host-only
+because libneo's metric is `class()`-dispatched and reads 3D splines, which
+cannot run under `!$acc routine seq`. `test/tests/test_cpp_canonical_device.f90`
+(built only with `SIMPLE_ENABLE_OPENACC=ON`) runs all three models for a batch of
+particles inside an `!$acc parallel loop` and checks the device result against
+the host step.
+
+`test/tests/test_cpp_canonical.f90` validates the analytic block against the
+regenerated python oracle. `test/tests/test_cpp_vmec.f90` runs the same
+integrator on `test/test_data/wout.nc` (a 2-field-period stellarator): the
+libneo metric satisfies `g g^-1 = I` to machine precision, `CP` energy stays
+bounded with no secular drift, and the big-step `CPP` orbit stays on a bounded
+radial band with radial bounce points, the guiding-center confinement signature.
+The stellarator is not axisymmetric, so the toroidal canonical momentum is not
+conserved and is not asserted; near the axis `s -> 0` the flux metric is singular
+and the central-difference gradients lose accuracy, so the test starts at
+mid-radius.
+
+The genuine 6D canonical CPP is wired into the production alpha-loss pipeline as
+`orbit_model = ORBIT_CPP6D` (5), through `COORD_VMEC`. `init_cpp` in `simple.f90`
+replicates the `init_sympl` sqrt(2) block (`mu` by factor 2,
+`ro0_bar = ro0/sqrt(2)`, `vpar_bar = vpar sqrt(2)`), reading `|B|` and `h_i` from
+the native VMEC field at the start, then seeds the state at `(s, theta, phi)`
+(s direct, no rho) with `MODEL_CPP_SYM`, `mass = 1`, `dt = dtaumin/sqrt(2)`, and
+`rho0 = ro0_bar`. Keeping `mass = 1` is what makes the wire well-posed: with the
+consistent `|h|^2 = 1` metric the kinetic term `(1/2m)(p-qcA)g^ij(p-qcA)` along
+the field reduces to `vpar_bar^2/2`, so the 6D Hamiltonian equals the GC one and
+the velocities stay `O(vpar_bar)` (physical-CGS `mass ~ 1e-24` would blow up
+`v^i = g^ij(...)/m` and wreck the Newton). The covariant momenta
+`p_theta = vpar h_theta + A_theta/ro0_bar`, `p_phi = vpar h_phi + A_phi/ro0_bar`
+match the GC seed; `p_s` carries the `g_si v^i` metric term, the genuine 6D
+start. `orbit_timestep_cpp_canonical` advances one `dtaumin/sqrt(2)` step and
+writes `z(1:5)` itself (`z(1) = s` for `COORD_VMEC`, `z(1) = rho^2` for
+`COORD_CHARTMAP`, angles direct, `z(4) = pabs`,
+`z(5) = vpar_bar/(pabs sqrt(2))`), so `times_lost`, `confined_fraction`, and the
+trajectory output read `z(1:5)` exactly as on the GC path. The macrostep in
+`simple_main` dispatches on `orbit_model`; the GC default and `ORBIT_PAULI` keep
+`to_standard_z_coordinates`, `ORBIT_CPP6D` routes around it. `ORBIT_CPP6D`
+requires a VMEC-backed canonical field (checked once in `main`, rejecting a
+standalone Boozer-chartmap input), needs `swcoll = .false.` and no wall, and
+error-stops in classification (`ntcut > 0` / `class_plot`); collisions, the wall
+path, and the classifier stencil are the documented follow-ups. The COORD_VMEC
+metric is attached once before the parallel region so per-thread `init_cpp` never
+races on the allocation.
+
+`test/tests/test_cpp6d_vs_gc.f90` drives the production `init_field` (BOOZER on
+the real `test/test_data/wout.nc`), then runs the production GC
+(`orbit_timestep_sympl`) and the production 6D wrapper (`COORD_VMEC`) from the
+same start over 2000 bare GC macrosteps. It asserts the chart consistency
+(`h_i g^ij h_j ~ 1`, the check the chartmap fails), `mu` held fixed and energy
+bounded, both orbits confined with overlapping radial (`s`) bands, and the
+`s >= 1` loss propagating through the wrapper. End to end, a 32-particle loss run
+on the same equilibrium gives a CPP6D confined fraction within `O(rho*)` of the
+GC one (`0.91` vs `0.97` at `trace_time = 1e-3 s`), the genuine-6D
+finite-Larmor difference over the GC.
+
+The full classical charged particle is wired the same way as
+`orbit_model = ORBIT_CP6D` (6), through `COORD_VMEC`. `init_cp` reuses the
+`init_cpp` sqrt(2) block and normalization (`mass = 1`, `qc = sqrt(2)/ro0`,
+`dt = dtaumin/sqrt(2)`) but seeds `MODEL_CP`, which resolves the gyration and
+drops the `mu|B|` term. So it needs the FULL velocity, not just the parallel
+piece: `v^i = vpar_bar h^i + vperp e_perp^i` with `vperp = sqrt(2 mu_bar |B|)`
+and `e_perp` a fixed-gyrophase metric-unit direction perpendicular to `h` (the
+raised radial direction with its `h`-parallel part removed). The seeded
+gyro-center sits an `O(rho*)` finite-Larmor offset off the GC start; that offset
+is the physics. `p_i = g_ij v^j + A_i/ro0_bar` as for `MODEL_CPP_SYM`. On the
+diagonal tokamak (`h_1 = 0`) the perpendicular direction reduces to the bare
+radial seed `[vperp, 0, 0]`, so the `COORD_TOK` oracle still reproduces bit for
+bit.
+
+Because the gyration is resolved, `ORBIT_CP6D` cannot run the bare GC macrostep;
+it must oversample the gyroperiod. The canonical cyclotron frequency is
+`Omega = qc |B| = |B|/ro0_bar`, so the gyroperiod in normalized tau is
+`2 pi ro0_bar/|B|`, while the step is `dtaumin/sqrt(2) = 2 pi rbig/(npoiper2
+sqrt(2))`. Steps per gyration is therefore `npoiper2 ro0/(rbig |B|)`, i.e. the
+resolution scales as `1/rho*`. On `test/test_data/wout.nc` (`ro0 = 2.70e5 cm`,
+`rbig = 1.02e3 cm`, `|B| = 5.60e4 G`, `rho* = ro0/(rbig |B|) = 4.7e-3`, gyroperiod
+`21.4 tau`) a sweep of `npoiper2` shows the energy error fall monotonically as
+the gyration is resolved: `max|dE/E0| = 0.31, 0.11, 0.033, 0.0087, 0.0022,
+0.00056` for `npoiper2 = 512, 1024, 2048, 4096, 8192, 16384` (2.4 to 77
+steps/gyration). It crosses below `1e-3` at `npoiper2 = 16384` (77
+steps/gyration), the required resolution there. A W7-X-class reactor case has a
+similar `rho* ~ 1/200`, so the same `npoiper2 ~ 1.6e4` order holds. At that
+resolution a single CP orbit conserves energy (`max|dE/E0| = 6e-4`), its
+gyro-center (running gyro-average of `s`) tracks the GC surface to `O(rho*)`, and
+the gyro-averaged magnetic moment shows no secular drift (the instantaneous
+`mu = vperp^2/(2|B|)` breathes at the gyrofrequency, `~8%`, which is not the
+invariant). `ORBIT_CP6D` shares `init_sympl` seeding,
+`orbit_timestep_cpp_canonical` (it dispatches on `cpp%model`), the loss
+write-back, and the `swcoll`/wall/classification guards with `ORBIT_CPP6D`;
+`test/tests/test_cp6d_vs_gc.f90` runs the sweep and the validation gates.
+
 ---
 
 ## 7. libneo Integration
@@ -885,6 +1103,10 @@ trajectory.
 | `src/orbit_symplectic_base.f90` | Integrator types and RK coefficients |
 | `src/orbit_symplectic.f90` | Symplectic methods |
 | `src/orbit_symplectic_quasi.f90` | Quasi-symplectic and RK45 |
+| `src/orbit_rk_core.f90` | Shared device LU and Newton shell |
+| `src/orbit_cpp_canonical.f90` | Curvilinear 6D canonical-midpoint integrator (cp/cpp_sym/cpp_var) |
+| `src/orbit_cpp_vmec_metric.f90` | VMEC metric/Christoffel + native field provider for the 6D integrator |
+| `src/orbit_cpp_chartmap_metric.f90` | Production Boozer/chartmap metric + field_can provider for the 6D integrator (ORBIT_CPP6D) |
 | `src/alpha_lifetime_sub.f90` | orbit_timestep_axis |
 
 ---