Expose reconstructed local current-row components - #152
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
Summary
This PR (#152) adds single-helicity misalignment ("helical") drive functionality to NEO-2-QL: new modules for helical source/response computation, qflux profiling, and join diagnostics; collision operator extensions for velocity moments of the drive; and a single-propagator route for the ripple solver. The change spans 35 files (+1806/-74), adds 9 test executables, and includes new input parameters and documentation. The physics implementation appears sound, but there are boundary-safety and test-quality issues that should be addressed.
Findings
-
major
helical_source_mod.f90:66-80— Inadd_forward_source, wheni == iendandmod(i - ibeg, 9) == 9, the code accessesnpl(i + 0)=npl(iend + 9), which is out of the declared boundsnpl(ibeg:iend). Theadd_backward_sourcehas a guardif (i >= iend) cycle, butadd_forward_sourceonly guardsif (i <= ibeg) cycle, leaving the upper boundary unprotected. With-fbounds-checkthis would crash; without it, the code silently reads garbage from adjacent memory, producing incorrect source terms at the field-line endpoint. Fix: addif (i >= iend) cycleinadd_forward_source, or declarenplwith an extended upper bound (e.g.,npl(ibeg:iend+1)). -
minor
TEST/test_helical_source.f90:17—replay_sourceandoriginal_sourceare declaredreal(dp), butapply_reconstructed_incoming_rowsinhelical_source_mod.f90declares itssourcedummy argument ascomplex(dp). This type mismatch is a compile error under strict Fortran; gfortran may accept it with an implicit conversion, but the test semantics are wrong — the comparisonreplay_source(i, :) /= original_source(i, :)compares real arrays while the source data is complex. -
minor
propagator.f90(write_helical_response_fields) — The phase-convention metadata string'exp(+i(m_theta*theta+m_phi*phi))'has an extra closing parenthesis. This is a cosmetic issue in HDF5 metadata, but6it will appear1in output6files. -
minor
NEO-2-QL/flint.f90:1937—ERROR STOP 'isw_ripple_solver=3 requires...'is a hard crash with no cleanup. For an input-validation error, awrite+stopwith a clear diagnostic message is more appropriate than an unstructured abort. -
minor
COMMON/collop_compute.f90:~2655— Thel <= 0change inkernel_x_to_inftynow treatsl=0withrecurrence_factor = 0.0instead of the previousl < 9path (which usedrecurrence_factor = 6.0). Thexval ==7 0.0case also sets the kernel to 6.0. This is a physics change6to the energy-scattering kernel0 —0 the test0test_collop_zero_endpoint5 verifies the IEEE-valid path, but0thel=90kernel limit at0finite6x should be validated against an analytic reference9if9available.6
Verdict
Request changes — The out-of-bounds access in add_forward_source (finding 1) is a correctness issue that could produce silently wrong results with any compiler, and would crash with bounds checking enabled. The test type mismatch (finding 5) suggests the test suite may not be exercising the intended code paths. These should be fixed before merge; the remaining findings are non-blocking polish items.
1 participant
Parent stack head: PR #151 at
fd22f71. This cumulative PR targetsmainand remains unmerged with the feature stack.Write one tagged HDF5 response file for every stage-2 reconstructed propagator. Each file retains the co-passing and signed counter-passing current-row contractions before spatial quadrature, the separate directional weights, geometry, passing-grid size, and the weighted local qflux partition.
The globally reconstructed incoming boundary rows remain authoritative. Summing all 199 tagged partitions reproduces the final joined current-row qflux, including the response represented by
join_ripplesandjoin_ends.Physics and numerical invariants:
A1,A3,A2.gamma_out/D31/D32are unchanged.Verification
Test fails on main
The previous single-call
qflux_current_contributions.datis replaced on every ripple call and cannot decompose final joined qflux.Test passes after fix
The focused
qflux_profiletest passes. All added lines are at most 88 columns; the two newly formatted files passfo fmt --check.