Validate the direct-EQDSK local field-line pitch - #91
Conversation
…d line PR #91 substitutes the local pitch hctrvr(2)/hctrvr(3) for the global q in the field-line label, but nothing tested that the substituted quantity is the pitch. The oracle is deliberately not the tabulated q: it is a field line integrated in cylindrical coordinates from field_eq alone, dR/dphi = R B_R/B_phi , dZ/dphi = R B_Z/B_phi , advanced until the geometric angle about the magnetic axis completes one circuit, giving the rotational transform as Delta_phi/(2 pi). That shares no code with do_magfie_eqdsk's coordinate transformation, so it is an independent behavioral oracle rather than a restatement of the patch. The two agree to 7e-4 at s = 0.15, 0.35, 0.60. The test also asserts the local pitch genuinely varies with theta. Without that an implementation returning the global q would satisfy the averaged comparison, and the test would not discriminate the change it exists to protect. Scope: this validates the substituted quantity, not the full canonical precession. The straight-field-line angle map also depends on s, so its gradient carries a radial term that the inp_swi=11 branch still omits; see #85. Writing this test surfaced a separate defect, now fixed in libneo: geoflux built both q and the s label from the GEQDSK's stored qpsi record, which on circ.eqdsk disagrees with the rotational transform of the fixture's own field by a factor of about 1.55. That is why test_eqdsk_axis could only afford a factor-of-two bound on the same comparison. This test needs the libneo fix to pass.
|
Update: added a test, and found the reason one could not be written before. The substituted quantity is now validatedThe change substitutes advanced until the geometric angle about the magnetic axis completes one circuit, The test also asserts the local pitch genuinely varies with theta — without that, Why the obvious test could not be written against
|
…ed qpsi (#404) A GEQDSK's stored `qpsi` record is not required to be consistent with its own `psirz` and `fpol`, and frequently is not — it may be a fit, sit on a different radial grid, or simply be stale. It should not be trusted for anything quantitative; q has to be computed from the field. `geoflux` trusted it for **two** things, so an inconsistent record corrupted more than the q profile: 1. the q returned by `geoflux_get_flux_profiles`, and 2. the toroidal flux, and hence **the `s` label itself** — `integrate_toroidal_flux` builds `s` by integrating `q dpsi_pol/(2 pi)`. So a bad `qpsi` means every direct-EQDSK quantity is evaluated at the wrong radius *and* labelled with the wrong q. ## How wrong, measured On `POTATO/test/golden_record_resonance/circ.eqdsk` the stored q disagrees with the rotational transform of the fixture's own field by a factor of about **1.55**: | s | q (stored, before) | rotational transform of the field | ratio | |---|---:|---:|---:| | 0.15 | 1.875 | 1.191 | 1.574 | | 0.35 | 2.292 | 1.537 | 1.491 | | 0.60 | 2.708 | 1.976 | 1.371 | This is also why `test_eqdsk_axis` could only afford a permissive factor-of-two bound when comparing `<B^phi/B^theta>` with q: a tighter bound would have failed on the fixture. ## The fix In axisymmetry with `B = F grad(phi) + grad(phi) x grad(psi)` and `psi` the poloidal flux **per radian**, `B_phi = F/R` and `B_pol = |grad(psi)|/R`, so a field line advances `dphi/dl_p = F/(R |grad(psi)|)` per unit poloidal arclength and one poloidal circuit gives ``` q(psi) = (F(psi)/(2 pi)) oint dl_p / (R |grad(psi)|). ``` Only `psirz` (through `psi_from_position`) and `fpol` are used, so the result is consistent with the field the rest of the chart evaluates. `locate_flux_surface`'s bisection is split out as `locate_by_normalized_flux` so contours can be traced from `psi` alone — necessary because `s` does not exist yet at the point q is needed. ## The two ends, which both need care - **Axis.** Surfaces too close to it collapse to a point and integrate to zero. The first *resolved* index is found rather than assumed to be 2: a zero left in place propagates into the toroidal flux and makes `s_grid` non-monotonic, which then corrupts every later interpolation. (This bit me during development — it produced `q ~ 1e12`.) - **Edge.** The outermost grid point sits on `sibry`, where the tracer has nothing to bracket against, and where a diverted equilibrium's q genuinely diverges. Left in, it produced a single `1.5e13` outlier that blew up `psi_tor_edge`. Both are extrapolated linearly in `psi` from the nearest resolved interior surfaces. ## Verification, against an independent oracle A field line integrated in cylindrical coordinates from `field_eq` alone, advanced until the geometric angle about the magnetic axis completes one circuit. It shares no code with the chart's coordinate transformation or with the q path. | s | q (computed, after) | field-line transform | ratio | |---|---:|---:|---:| | 0.15 | 1.19391 | 1.19077 | 1.0026 | | 0.35 | 1.53716 | 1.53669 | 1.0003 | | 0.60 | 1.97419 | 1.97610 | 0.9990 | From a factor 1.58 discrepancy to 0.03–0.3%. Note the ratio lands at `1.000` and not at `2 pi` or `1/(2 pi)`, which is the direct check that the flux-per-radian convention in the formula above is the right one — the usual place this calculation goes wrong. The matching regression test lives in NEO-RT ([itpplasma/NEO-RT#91](itpplasma/NEO-RT#91), commit `71c39fa`) because that is where the consuming chart is; it fails without this change. All 14 NEO-RT tests pass with it. ## Risk This changes q, and therefore the `s` label, for **every** direct-GEQDSK consumer. That is the point — the previous values were wrong wherever a file's `qpsi` was inconsistent — but it will move existing results, and anything that was tuned against the old labelling will shift. Worth a look from whoever owns the GEQDSK consumers.
![slopqueue[bot]](https://avatars.githubusercontent.com/u/48486978?s=60&v=4){kind=small}
There was a problem hiding this comment.
Review verdict: Approve
Summary: PR #91 fixes the canonical toroidal precession (Om_tB) in the direct-EQDSK (inp_swi=11) path of src/orbit.f90 to use the local field-line pitch B^phi/B^theta (from hctrvr(2)/hctrvr(2) (global, tabulated) — because in geoflux coordinates theta is the geometric poloidal angle, not a straight-field-line angle, so alpha = phi - q*theta is not a field-line label. A new test (test_eqdsk_local_pitch.f90`) validates the theta-averaged local pitch against an independent cylindrical field-line integration. The change is small, well-scoped, and well-documented.
Findings:
- minor
src/orbit.f90:447—q_local = hctrvr(3)has no guard against a vanishing denominator, whereas the test (test_eqdsk_local_pitch.f90:90) explicitly guardsabs(hcon(3)) <= tiny. For the tested fixture and any real tokamak ats > 0,B^thetais non-zero so this never triggers, but if a particle orbit ever reaches a near-axis region it would emit NaNs into the orbit integrator. Optional defensive guard (or at least a documented precondition) would harden it; not a blocker.
Verdict: Approve — the physics fix is correct (the local pitch is the actual dphi/dtheta = B^phi/B^theta), the test is a genuine independent oracle (cylindrical field_eq shares no code with do_magfie_eqdsk's metric path) and is non-vacuous (it asserts the local pitch varies along theta, so an implementation returning the global q would fail), and the honestly-scoped limitation (the radial gradient term of #85 is still omitted) is clearly documented in both code comments.
krystophny
left a comment
krystophny
left a comment
There was a problem hiding this comment.
Reviewed latest head 85d6ef8 against updated #90 after merging current main. The change uses the local direct-EQDSK field-line pitch B^phi/B^theta where geometric theta is not a straight-field-line angle. The new regression is a genuine independent oracle: it integrates the cylindrical field directly, requires non-constant local pitch, and checks the theta average. The full integrated fo pipeline passes.
The prior minor hardening note still applies: production divides by hctrvr(2) without the near-zero guard used by the test. That is not a blocker for the validated off-axis tokamak domain, but the precondition should eventually be explicit if near-axis orbits are admitted. No blocking finding.
krystophny
left a comment
krystophny
left a comment
There was a problem hiding this comment.
Reviewed latest cee794ec relative to updated #90. The local field-line pitch remains validated against independent cylindrical field-line integration, and the prior denominator hardening note is fixed: production now fails closed before dividing when direct-EQDSK B^theta vanishes. Full integrated fo is green. No blocking finding.
…d line PR #91 substitutes the local pitch hctrvr(2)/hctrvr(3) for the global q in the field-line label, but nothing tested that the substituted quantity is the pitch. The oracle is deliberately not the tabulated q: it is a field line integrated in cylindrical coordinates from field_eq alone, dR/dphi = R B_R/B_phi , dZ/dphi = R B_Z/B_phi , advanced until the geometric angle about the magnetic axis completes one circuit, giving the rotational transform as Delta_phi/(2 pi). That shares no code with do_magfie_eqdsk's coordinate transformation, so it is an independent behavioral oracle rather than a restatement of the patch. The two agree to 7e-4 at s = 0.15, 0.35, 0.60. The test also asserts the local pitch genuinely varies with theta. Without that an implementation returning the global q would satisfy the averaged comparison, and the test would not discriminate the change it exists to protect. Scope: this validates the substituted quantity, not the full canonical precession. The straight-field-line angle map also depends on s, so its gradient carries a radial term that the inp_swi=11 branch still omits; see #85. Writing this test surfaced a separate defect, now fixed in libneo: geoflux built both q and the s label from the GEQDSK's stored qpsi record, which on circ.eqdsk disagrees with the rotational transform of the fixture's own field by a factor of about 1.55. That is why test_eqdsk_axis could only afford a factor-of-two bound on the same comparison. This test needs the libneo fix to pass.
d97892a to
77f73f3
Compare
cee794e to
4454d14
Compare
krystophny
changed the title
1 participant
Summary
The direct core already projects drift with the local field-line pitch. This PR now contains the independent validation and fail-closed guard only:
B^phi/B^thetaB^thetabefore forming the local pitchThe superseded duplicate implementation commit was dropped during the stack rebuild.
Validation
fos=0.15,0.35,0.60Stack: #90 -> this PR.