Fix covariant B_theta and add psimax stop in EQDSK chartmap converter#399
Merged
Conversation
The converter wrote B_theta = G/q, reasoning "iota = I/G". That confuses the covariant with the contravariant component ratio: iota is B^theta/B^phi, while the covariant I(s) measures the enclosed toroidal current (2*pi*I = mu0*I_tor). For a tokamak I is ~40x smaller than G/q, and SIMPLE consumes B_theta as the covariant surface function (h_theta = B_theta/Bmod), so the fabricated value inflated the parallel field-line length factor (G + iota*I) and slowed every traced orbit's parallel dynamics by O(1/q^2): 16 % too-low bounce frequencies on the circular benchmark equilibrium (q = 2.26 at s_tor = 0.3). Compute I(s) instead as the signed circulation of B_pol along the flux surface, with B_pol = |grad psi|/R taken from the g-file psi map itself (sqrt(B^2 - B_tor^2) from the Boozer splines cancels catastrophically on inner surfaces). Verified against POTATO tracing the same EQDSK: bounce and toroidal frequencies of trapped and passing 5 keV deuterons now agree within 0.7 % (previously 16 % off). Also expose psimax (poloidal boundary flux, CGS) through convert_eqdsk_to_chartmap. Without it the flux-surface scan locates the boundary purely by field-line integration, which for synthetic equilibria without an X-point (psi rising to the box edge) runs to nsurfmax through hours of vacuum-region surfaces and labels the wrong surface s = 1. Tests: Ampere check of I at the outermost surface against the g-file Ip, |I/G| ordering, and psimax .inp pass-through.
krystophny
added a commit
to itpplasma/benchmark-simple-potato
that referenced
this pull request
Jul 16, 2026
Replace the R0 = 1.6 m circular equilibrium by a reactor-size case (R0 = 6.2 m, a = 2.0 m, B0 = 5.3 T) so that both the 5 keV deuterium benchmark particle and SIMPLE's default 3.5 MeV alpha are confined at mid-radius. gen_circular_eqdsk.py writes the EQDSK with a poloidal flux realizing q(r) = 1.5 + 2.5 (r/a)^2 exactly on the concentric-circle model. Add two conversion tools usable for any axisymmetric g-file: - tools/eqdsk_to_simple_chartmap.py: libneo EQDSK-to-Boozer converter plus SIMPLE regrid. Stops the converter's flux-surface scan at the g-file boundary flux (needed for equilibria without an X-point) and requires libneo with the corrected covariant B_theta (itpplasma/libneo#399, which fixes 16 % too-slow parallel dynamics from the earlier I = G/q). - tools/eqdsk_to_vmec.py: VMEC++ fixed-boundary re-solve with pressure and q profile derived from the psi map; produces rung0/wout_circ.nc for SIMPLE's native VMEC input path. Cross-code validation on the committed files (5 keV deuteron, s_tor = 0.3, canonical frequency API vs POTATO on the same EQDSK): omega_b and omega_phi agree within 0.7 % for trapped (xi = 0.3) and passing (xi = 0.8) orbits; the VMEC variant agrees to 2.5 % in omega_b (force-balanced re-solve, Shafranov shift). Checksums and the numbers are recorded in rung0/README.md. Drop the superseded small-machine chart maps, the nzeta expansion helper (folded into fix_chartmap_for_simple.py), and the POTATO-generated eqmagprofs.dat run product.
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.
Two defects in
eqdsk_to_boozer_chartmap(follow-up to #346), found while regenerating the reactor-size circular benchmark fixture for itpplasma/benchmark-simple-potato:1.
B_theta = G/qis not the covariant poloidal component. The inline comment reasonediota = I/G => I = G/q, but iota is the ratio of the contravariant components; the covariantI(s)measures the enclosed toroidal current (2*pi*I = mu0*I_tor). For a tokamak the trueIis ~40x smaller thanG/q. SIMPLE consumesB_thetaas the covariant surface function (h_theta = B_theta/Bmod), so the fabricated value inflated the field-line length factor(G + iota*I)and made every traced orbit's parallel dynamics too slow byO(1/q^2)— 16 % too-low bounce frequencies atq = 2.26. Now computed as the signed circulation ofB_pol = |grad psi|/Rfrom the g-file psi map (usingsqrt(B^2 - B_tor^2)instead would cancel catastrophically on inner surfaces).Cross-code check on the circular benchmark EQDSK (5 keV deuteron,
s_tor = 0.3), SIMPLE canonical frequency API on the converted chartmap vs POTATO adaptive-RK tracing of the same g-file:(before the fix: omega_b was 1.0796e4 / 3.1317e4, i.e. 16 % low)
2.
psimaxwas hardwired to the disabled placeholder 1e10. For real diverted g-files the field-line scan leaves the box just past the X-point, but synthetic equilibria without an X-point have psi rising smoothly to the box edge: the scan then grinds through thousands of vacuum-region surfaces (hours instead of seconds, sinceflint_for_Boozeralso ignores itsnstepargument) and labels the last box-confined surfaces = 1instead of the LCFS.convert_eqdsk_to_chartmapnow acceptspsimax(CGS, absolute) and passes it intoefit_to_boozer.inp; default behaviour is unchanged.Tests: Ampere check of
Iat the outermost stored surface against the g-fileIp,|I/G|ordering, andpsimax.inp pass-through. Fullpython/tests/test_eqdsk_to_boozer_chartmap.py+test_eqdsk_base.py: 13 passed.