port/c: full 9-case gate -- honest surface-dependent agreement - #165
Draft
krystophny wants to merge 17 commits into
Draft
port/c: full 9-case gate -- honest surface-dependent agreement#165krystophny wants to merge 17 commits into
krystophny wants to merge 17 commits into
Conversation
Literal C port of itpplasma spline (spline_coeff/spline_val_0), calling the same LAPACK dptsv so coefficients are bit-identical. Cross-checked against the real Fortran spline module: 163 checks (coeffs + value/1st/2nd derivative over 41 points), 0 failures at rtol 1e-12. Baseline: clean rebuild 16.1s, ctest 0.88s (3/3), neo_rt.x run 0.16s. Golden gate rtol=1e-8; golden-path deps reduce to spline + DVODE (no BLAS/LAPACK/NetCDF beyond dptsv inside spline).
The golden gate at rtol=1e-8 only holds single-threaded: under OMP_NUM_THREADS=1 the Fortran build reproduces its own golden.h5 bit-exact (9/9, max abs diff 0.0). Multithreaded it diverges at ~1e-4 (run-configuration dependent), which is why golden_record is excluded from ctest. All port branches gate single-threaded against this reference via run_gate.sh. golden.h5 stays gitignored and persists across branch checkouts.
Literal C port of do_magfie_mod: Boozer reader (inp_swi=8/9), spline-based profile/mode interpolation, and field evaluation. Cross-checked against the real Fortran do_magfie_mod over 33 theta points x 9 quantities (bmod, sqrtg, grad-B, hcovar, hctrvr): 297 checks, 0 failures at rtol 1e-12. Includes do_magfie_pert amplitude port (used next for the magfie.out bn column).
Measured: GSL msadams matches DVODE to 5.6e-9 rel at rtol=1e-9/atol=1e-10 on a smooth ODE -- under the 1e-8 gate but <2x margin. Orbit root-finding events are the residual risk; fallback is porting DVODE's Adams path. GSL 2.8 and SUNDIALS CVODE 7 both available.
Literal C port of neort_profiles (plasma.in/profile.in readers, plasma+rotation splines, thermodynamic forces) and collis_alp (loacol_nbi/coleff/onseff). Cross-checked against the real Fortran modules: 21 scalars (vth, M_t, Om_tE, A1, A2, ni/Ti/Te, qi, mi, efcolf/velrat/enrat) at rtol 1e-12, 0 failures. collis keeps its local ev=1.6022e-12 distinct from util ev=1.602176e-12, faithful to the Fortran. loacol dchichi/slowrate outputs are unassigned in Fortran (caller discards); zeroed here.
bounce_integral uses DVODE root-finding (nevents=2). GSL gsl_odeiv2 has no event location; SUNDIALS CVODE (CVodeRootInit, CV_ADAMS) is the faithful match.
SUNDIALS CVODE (CV_ADAMS + SUNNonlinSol_FixedPoint) reproduces DVODE to machine precision (9e-15 rel) at the code's tolerances, vs GSL's 5.6e-9 -- they share the VODE Adams lineage. CVODE also provides CVodeRootInit for the bounce-integral events. Decision: use SUNDIALS CVODE. Integrator removed as a project risk.
Orbit bounce integration ported to SUNDIALS CVODE (CV_ADAMS + fixed-point + CVodeRootInit) mirroring DVODE Adams + g_fcn=bounceroots. magfie flux-surface average (eps, B0, Bmin, Bmax, etatp, etadt, th0) ported. Cross-checked vs the real Fortran DVODE orbit on 3 real orbits (trapped + passing), bounce_time + taub + bounce averages: 18 checks, all within the 1e-8 gate. Worst case 8.07e-9 on a deep-trapped vpar average -- event-location (turning-point root) sensitivity, vs 9e-15 on a smooth proxy ODE. Integrator substitution validated on real physics.
Ported neort_freq (canonical-frequency splines, Om_th/Om_tB/Om_ph, d_Om_ds) plus orbit bounce_fast/orbit_get_s. Cross-check vs Fortran neort_freq: in-range frequency VALUES agree to 1.48e-8 -- the propagated CVODE-vs-DVODE orbit difference. Near-tpb extrapolation amplifies it to ~2e-2 (values) / 15% (derivatives) over a tiny eta sliver. Key finding: the 1e-8 golden gate effectively requires bit-reproducing DVODE. Faithful ports with CVODE realistically reach ~1e-7. Decision needed on gate tolerance vs porting DVODE exactly. See PROGRESS.md.
Ported neort_resonance (driftorbit_coarse/nroot/root): resonance-line bisection on the freq layer. Added do_mph (toroidal mode) sourced from the perturbation file (field_pert_mph = nint(nfp*n)). All 9 C modules compile clean under -Wall -Wextra. Resonance is pure bisection on freq; validated end-to-end via the executable gate at rtol 1e-7.
Ported neort_transport (compute_transport_integral midpoint rule, D11/D12/Tphi integrands, timestep_transport with perturbation Hamiltonian) and neort_nonlin (nonlinear_attenuation: exact 1.0 for golden nonlin=false; omega_prime ported; nonlin=true attenuation_factor path errors loudly, outside golden gate). Generalized orbit bounce_fast_ext(custom RHS, istate) + exposed poloidal_velocity. Mirrors the Fortran 'if(nroots==0) continue' no-op gotcha. All compile clean.
Driver (namelist parser, neort_compute_at_s + compute_transport orchestration, check_magfie, 5 output writers). neo_rt_c.x builds and runs the golden path. Case 0p500 vs golden: magfie bit-perfect (0.0), dominant resonance mth=5 matches to 7 digits. Discrepancy concentrated in near-threshold trapped resonances (mth=1..3 trapped ~7-13%) where freq log-extrapolation from near-separatrix bounce points amplifies CVODE-vs-DVODE divergence. Full-transport agreement ~6e-4 on totals -- coarser than the 1e-7 in-range estimate; near-separatrix resonances are the limiter. Run ~4.6s vs Fortran ~6.5s single-threaded.
Reusable compare_golden.py (any port exe, configurable rtol). Full 9-case run corrects the earlier 0p500-only estimate: cross-language agreement is surface dependent -- dominant D11 within ~1e-3 at low s but 2.5-5% at high s (0p700-0p900), sub-dominant harmonics far worse. 0p500 was the best case, not representative. Higher-s dominant resonances sit in the near-separatrix region where CVODE and DVODE diverge most. Field/magfie output stays bit-exact.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This draft PR publishes the existing committed state of
port/cso the local branch/worktree can be retired after review.Checks were not run in this cleanup operation; the branch-specific CI should provide validation.