Preserve double precision in VACUUM observer-grid construction

[krystophny]

VACUUM constructs its Cartesian observer grid through single-precision temporary arrays and single-precision grid-spacing expressions. This loses coordinate precision even though the caller and stored observer coordinates use double precision. The patch changes the two temporary-array declarations to REAL(r8) and the two spacing constants to 1.0_r8.

Thanks to Nik and Jake for explaining the near-surface quadrature limitation. This four-line patch preserves the existing snapping, surface traces, evaluator, public API, and output formats. Grid resampling belongs in downstream postprocessing, where snapped or unsupported samples can be masked without accepting inaccurate near-surface evaluations. The separate field-content limitation in issue 171 remains open.

Downstream native-cell resampling is available in rmp_torque, with explicit validity masks and separate inside/outside interpolation. TC24 pins the same feature on top of its existing converter lineage.

Validation against upstream develop at 5be646f2d67f0d6002c531ff2b855e701103bbec: the independent weighted-endpoint grid oracle fails on the base at 4.30e-8 relative error and passes this patch at 1.48e-16. Independent inside/outside surface-trace checks pass on both versions, with identical snapped coordinates, flags, and printed fields. The native GPEC/DCON build and both registered focused CTests pass. The regression harness is retained outside this PR; the upstream diff contains only these four precision changes. Gremlin remains blocked before capture by the existing declared-Git-dependency tool defect; these results come from native CMake/CTest.

Chris&AI

[logan-nc]

@claude review this

[claude]

Claude finished @logan-nc's task in 42s —— View job

---

Review of PR #294: Cartesian vacuum query coordinates

I reviewed vacuum/vacuum_ma.f and gpec/gpout.f by reading the diff. I did not build the code or run the tests, so everything below is from the source. I read the tests/ harness only at a high level, and I skipped the CMake and Python fixture files.

Summary: The opt-in design is sound.

• Absent or false preserve_query_coordinates leaves the loops path and the snapping branch unchanged.
• The exact path feeds double-precision query_r and query_z straight into xloop and zloop.
• It validates shapes, the grid dimensions and nloopr == 0 before use.
• Unavailable queries get igdl = -2, complex NaN outputs, and are skipped in chi through query_valid.
• All four Cartesian mscfld calls in gpout_xbrzphi forward the coordinates. The non-Cartesian call stays legacy.

Issues

1. Lines over 72 columns in a fixed-format file (CLAUDE.md requires ≤72).

   • vacuum/vacuum_ma.f:744 is code, and it is the one that matters:
     $ ERROR STOP 'Exact pickup requires a nonzero derivative step'
     It is past column 72, so a strict fixed-form build truncates the string literal and fails to compile. It may compile today only if the build passes -ffixed-line-length-none or similar. Please split it with // concatenation or shorten the message.
   • vacuum_ma.f:717 and :997 are new comments over the limit. They are harmless but should be wrapped.
   • Lines 571 and 813 were already over the limit in the base commit, so they are not caused by this PR.

2. query_valid is not initialised on the legacy path in the way chi expects.

   • query_valid is set to .true. for all ndimlp entries and then passed to chi in both the legacy and exact paths.
   • This is correct, but it means legacy behavior depends on that blanket initialisation. A short comment would guard against someone removing it.

3. Non-positive R test. IF (xloop(i) <= 0.0_r8) masks queries with R ≤ 0. Is that intended for every geometry? It is reasonable for a tokamak grid. The README or a comment should say it is a deliberate part of the contract.

4. Shared guard scope.

   • check_cartesian_postprocessors runs on btr/btz/btp after assembly, which is the right place to catch NaNs before gpeq_rzpdiv or Chebyshev processing.
   • With both flags off, NaNs propagate to the raw export. That is documented in the comment, but downstream NetCDF/ASCII writers will now contain NaNs. Please confirm no reader of those outputs assumes finite values.
   • The guard sits inside the same branch as the first mscfld pair. Please confirm the other gpout_xbrzphi paths that call mscfld (around lines 4749-4770) have equivalent protection, or are not followed by divzero or Chebyshev.

5. Clearance semantics.

   • clearance = 2*max(|delx|,|delz|) is measured against the discretized xpla/zpla and xwal/zwal contours.
   • Queries inside that band become unavailable, not just those on the contour. This is the intended "unresolved stencil" behavior.
   • It does change which grid points are valid compared with the old snapping, which tagged them -1. Please note this in the docs or release notes, since users will see new NaN holes near the plasma boundary and wall.

Minor

• ERROR STOP inside library routines aborts the whole run. That is acceptable for contract violations, but it is harsh for the zero derivative step case.
• bxr(i) = (chir(3,i)-chir(4,i))/(2.0*delx) still mixes a default-real 2.0 with double-precision operands. This is pre-existing, but you could change it to 2.0_r8 while you are in the area.
• The PR description says 10 files and +617/-13. That matches the commit.

Not verified

I did not review the test logic in test_vacuum_query.f90 or test_cartesian_postprocessors.f90 in depth. I did not check CMake portability or run the fo pipeline.

Recommendation: Fix item 1 (the line-length problem at vacuum_ma.f:744) before merging. Items 3-5 are worth a short reply or a doc note.

[logan-nc]

@krystophny the llm says your implementation is sound, but as a human I see a lot of new code spent changing what is an intended behavior with no clear use case motivating so much change. The grid is intentionally snapped because the Fortran vacuum code gives incorrect numbers for near-surface evaluation points. This is something @jhalpern30 is improving in his Julia implantation. The GPEC behavior has been around forever and is the official io... the primary approach to using gpec should be to use this irregular grid in your post processing. You can interpolation regular grid from there if you wish. If you have a use case where you need the incorrect singularity-caused garbage at points near the control surface please explain the workflow a bit to us in human terms. In my experience, the post-processing is worse this way.

If it's just the real4 rounding of unsnapped points away from the controlsurface that's the issue, i think this could be a much more simple fix to move those to real8 without requiring a new opt-in flag. I fully support that.

The new test suite seems like it would be a good direction to go if we planned a long term modernization of the Fortran code but is perhaps a bit odd to add here in a limited capacity just before the final Fortran release (1.6) is expected. I'd be curious what @matt-pharr and @parkjk think about the pros/cons of it.

[jhalpern30]

Agreed with @logan-nc regarding the point-snapping - the vacuum code can't properly evaluate the vacuum fields just outside but not on the surface, it can only do exactly on the surface or sufficiently far from it, hence the snapping behavior. It seems that this PR would allow (but not fix) evaluation much closer than currently allowed, which I think @logan-nc has shown can already be rather finicky in its current state.

And agreed that interpolating these values on whatever grid you want in postprocessing is probably the way to go. At some point I'll be implementing the near-singular quadrature in Julia, so if you have a motivating example that needs precise fields very close to the surface, I'd be happy to do it sooner rather than later.

[krystophny]

Thanks for checking this. I reduced the code changes to the minimum and do the rest in post-processing now.

[matt-pharr]

@logan-nc I agree with your comments about tests + additions. Let's keep it minimal. This PR is good as-revised.