Reject numerical ties in field-line return selection - #153
Conversation
isw_qflux_NA=1 with in_file_pert='none' assembles an identically zero non-axisymmetric source. That zero vector is the Arnoldi seed, so its normalization divides by zero (0/0 -> NaN) and zhseqr then fails to converge at every step, spinning without qflux output. Reject the configuration in check() with a clear message instead.
Add an optional right-hand side to the non-axisymmetric solve: the drive of a single-helicity perturbation given by a radial corrugation field (delta B^s/B_0^phi)_mn and an electrostatic potential harmonic Phi_mn. The vpar corrugation piece enters the A1/A2 force channels with velocity weights a3m and x^2 (sigma-odd), the ExB piece with weights x^-1 and a1m (sigma-even); the two new velocity moments are computed in compute_sources. Input via new &ntv_input variables (isw_hel_drive, m_theta_hel, hel_brad_re/im, hel_phim_re/im), default off. Band profiles, moments and dimensional prefactors follow a CAS-checked derivation (see issue #118). Off by default: golden-record QL summary.h5 is reproduced to 1e-15 against a main build.
Review findings on the drive commit: fold the x^(3k-4) weight into a single power so the x=0 quadrature node stays finite; skip the moment quadrature on the precomputed-matrix-element path and reject isw_hel_drive at run time when the moments are unavailable; reject the unvalidated relativistic branch, reconstruction mode and unset species data; deallocate asource_hel in collop_unload.
The zero-source guard added upstream rejects isw_qflux_NA=1 with in_file_pert='none'. With the misalignment drive, that configuration is valid: the drive supplies the non-axisymmetric source. Relax the guard so it fires only when the drive is off or all its amplitudes are zero.
The opposite-helicity drive (m,n) -> (m,-n) enters add_helical_source as the complex conjugate of the band profile (phase convention exp(+i(m theta + n phi))). The assembly stencil is a real-coefficient linear map, so it must conjugate the assembled source pointwise; a sign or phase error in the stencil breaks this. The existing checks use a profile whose real and imaginary parts vary independently, so they do not exercise conjugate symmetry. Add a check that conjg(profile) yields conjg(source) to 1e-14, plus a vacuity guard that the source carries an imaginary part. Verification: - ctest helical_source_test: Passed (0.07 s) - ctest -LE external-data: 18/18 passed - red-before: injecting a complex coefficient (1,0.1) into one stencil term makes the new check error stop with 'conjugate profile did not conjugate the source pointwise' (exit 1)
Decompose the phi-integrated current channel qflux(2,k) = sum_col flux_vector(2,col)*source_vector(col,k) into its per-field-line-point contributions in the mode-1 local (axi) solver, which HELCORE uses. The new pointwise_current_mod reuses the already-verified flux_vector current row: its counter-passing block enters with the explicit minus sign (parallel current is odd under sigma -> -sigma) and the step_factor phi-quadrature weight is baked in, so the decomposition introduces no new sign, weight, or normalization. flux_vector's columns partition into disjoint per-istep blocks, so the block dot products sum back to qflux(2,k) exactly. The solver writes current_profile.dat (phi_mfl, per-drive contributions, and current densities) guarded by isw_hel_drive, so runs without the helical drive are bitwise unchanged. A hard self-check stops the run if the field-line sum does not reproduce qflux(2,k). Verification: - 17/17 TEST ctests pass; new pointwise_current_test added. - Red-before: breaking the block partition (nblk without the npassing+1) gives "ERROR STOP FAIL: per-point contributions do not sum to the qflux current channel"; restored -> Passed. - Local HELCORE simple_l1 run, all three drive cases (nminus_phi, nminus_phi_im, nplus_phi): gamma_current/D31/D32 reproduce the accepted aCluster response (-7.338706/-20.944508, -2.711060/-3.138948, -2.369674/-1.825414) to ~1e-11; each current_profile.dat reconstructs qflux(2,k) to ~1e-15 relative over 961 finite points.
![slopqueue[bot]](https://avatars.githubusercontent.com/u/48486978?s=60&v=4){kind=small}
There was a problem hiding this comment.
Review verdict: Request changes
PR #153 Review
Summary: This PR adds a single-helicity misalignment (helical) drive source to the QL NTV solver, introduces a source-only non-axisymmetric solve path (in_file_pert="none" with prescribed m_phi_input), refactors ripple-solver routing and join diagnostics into small modules, fixes a numerical singularity at the x=0 quadrature endpoint in the collision operator, corrects spline-shift loop bounds, adds a relative tie tolerance for field-line return distances, updates D_ij metadata units (cm² G/s for current channels), and adds 14 new unit tests plus CI test labels. The physics is new but self-contained and guarded by runtime switches defaulting to no-op for existing runs.
Findings:
-
[major] COMMON/collop_compute.f90 — The
am_helintegrand constructs weights proportional tox**(3*k + alpha)for the helical drive source. Whenalpha < 0(which the Maxwellian/Laguerre setup permits for the non-relativistic path used byisw_hel_drive), thex=0endpoint ofintegrate(am_hel, 0d0, ...)diverges, and unlike the existingkernel_x_to_inftypath (which special-cases x=0),am_helhas no x=0 guard. This NaN can propagate intoasource_heland then into the RHS of the non-axisymmetric solve. Either add analpha >= 0assertion on the hel-drive path, or special-case x=0 inam_helas is done forkernel_x_to_infty. -
[minor] NEO-2-QL/neo2.f90 —
has_perturbation_file()now gates the symlink creation forin_file_pert. This is correct, but the function returns.false.whenin_file_pert == "none", which is now the default. Existing users who relied on the old default (noin_file_pertin their namelist) will now silently skip the perturbation file symlink. This is intentional for the new source-only path but should be noted in release notes. -
[minor] NEO-2-QL/neo2.f90 — The new
check()guard errors out whenisw_qflux_NA == 1and no perturbation file and no helical drive source are present. This is good defensive coding. However, the error message says "isw_qflux_NA=1 with in_file_pert=none provides no perturbation field" which is accurate, but the testntv_driven_source_guard_testonly covers the "driven" path. Consider adding a test for the "no source at all" rejection path to lock in the error message.
Verdict: Request changes — The am_hel x=0 endpoint (finding 1) needs an explicit guard to prevent NaN propagation into the non-axisymmetric solve; the existing kernel_x_to_infty precedent shows this is a known numerical issue that should be handled consistently.
1 participant
Summary
The change leaves field integration, return-distance evaluation, propagator construction, ABI, and memory layout unchanged. It only changes the comparison used to select the terminating return.
Verification
Failing before, with fixed geometry and increasing steps per period:
The 3840-step result selected a numerical near-tie that was 3.5E-9 smaller than the period-1 distance.
Passing after:
The formatter reports the existing legacy whole-file style warning for
COMMON/mag_interface.f90; no broad reformat is included.