Resistive singular coupling

[matt-pharr]

Changes related to adding scientifically robust singular coupling to the RDCON -> GPEC workflow. WIP.

[matt-pharr]

This is what was accomplished by adding the resolution option for the inner layer (was previously limited to 3 points):

The plot on the bottom has 60 points in each resonant and extension cell instead of 3. This does not affect asymptotic matching, but does make plots of xi/Jbgradpsi look better.

[matt-pharr]

A scan of resistivities using a LAR equilibrium with only a q=2 surface shows interesting behaviors in xi. The direction of the interchange switches at eta=8e-8, or in this case for lundquist numbers lower than about 10^6.

Scan of higher lundquist numbers should be done.

[matt-pharr]

Choosing sing_spot based on the first plot posted in this pull request (which puts it just outside of the resonant layer using default rdcon dx/layer width values) metrics in a single-surface LAR case are relatively well-behaved.

B_pen is far more dependent on the resistivity / lundquist number in this case.

[matt-pharr]

@StuartBenjamin you are one of the few who knows the guts of RMATCH... can you take a very brief look at these comments (in f27f665) and let me know if there should be any changes/if you know what some of the 'unknown's are?

[StuartBenjamin]

Remembered a few of the extra terms!

[matt-pharr]

@logan-nc The problems wound up being that the way rdcon was regridding was cutting off at a different psilim than I expected, resulting in the poorly conditioned solution set that we know occurs at some psi cutoff values. I haven't run it yet, but I think I should be able to get it to work nicely with Stuart's regridding contributions. Hopefully I'm not dooming it by speaking too soon!

[matt-pharr]

@claude Please briefly summarize this PR.

[claude]

Claude finished @matt-pharr's task in 1m 9s —— View job

---

Summary of PR #210: Resistive Singular Coupling

This PR enhances GPEC's resistive MHD analysis capabilities by implementing scientifically robust singular coupling between RDCON (resistive stability) and GPEC (perturbed equilibrium) codes. The work enables accurate 3D field penetration calculations at rational surfaces where resistive effects dominate.

Key Improvements

1. Inner Layer Resolution Enhancement

• Previously limited to 3 grid points in resonant layers, now supports configurable resolution (e.g., 60 points)
• Improves visualization quality of ξ/J·∇ψ profiles near singular surfaces
• Physical resolution maintained consistently through width-aware grid scaling

2. Resistive Layer Coupling

• Fixed critical bug where RMATCH output mismatched psi points and solution points, causing spurious zeros in spline interpolation
• Now properly excludes only the singular point at rational surfaces (matching RMATCH behavior)
• Enables intra-inner-layer scanning of resonant quantities

3. Field Penetration Physics

• Demonstrated actual 3D field penetration at flux surfaces (previously broken)
• Resistivity scans show physically meaningful behavior: interchange direction switches at η≈8×10⁻⁸ (Lundquist number ~10⁶)
• B_pen (penetrated field) shows expected strong dependence on resistivity/Lundquist number

4. Singular Spot Positioning

• Implemented logic for choosing sing_spot (interpolation point outside resonant layer)
• Well-behaved metrics achieved when positioned just outside default RDCON layer widths

5. Code Infrastructure

• Refactored RDCON/RMATCH I/O to reduce write overhead and improve data handling
• Added NetCDF outputs for resistive quantities, Lundquist number, and sweet spot metrics
• Improved RMATCH parameter documentation (with help from @StuartBenjamin)
• Separated equilibrium namelist loading for cleaner code organization
• Added XDRAW visualization files for RPEC (resistive perturbed equilibrium) solutions

6. Numerical Robustness

• Resolved regridding issues with psihigh that were causing poorly conditioned solution sets
• Improved deltac numerics and cubic Hermite finite element method parameters
• Fixed real→complex eigenvalue conversion from RMATCH

Files Changed

100 files modified across RDCON, RMATCH, GPEC, EQUIL, and example input files (+2613/-1894 lines)

Status

Work-in-progress with ongoing refinements to resistive layer width logic and interpolation scaling. Recent commits focused on finalizing parameters for resistive metrics publication.

--- | Branch: RDCON_singcoup

[logan-nc]

@matt-pharr can you confirm the truncation issue is fixed and this is ready for full review? Last message from you to me ended on a cliff hanger, "Hopefully I'm not dooming it by speaking too soon!"

[logan-nc]

@claude review this PR

[claude]

Claude finished @logan-nc's task in 6m 29s —— View job

---

PR #210 Review: Resistive Singular Coupling

Tasks:

• Gather PR context
• Review core RDCON/RMATCH changes (rdcon/, gpec/rdcon.f)
• Review GPEC output and coupling changes (gpec/gpout.f, gpec/gpec.f)
• Review equilibrium changes (equil/)
• Review input file changes
• Review numerical robustness and physics correctness
• Post final review

---

This is a substantial WIP PR implementing resistive singular coupling. The physics approach (using RMATCH eigenvalues and layer widths to set per-surface sing_spots, reading globalsol.bin with correct format, interp_np_res for inner-layer resolution) is sound. However, there are several bugs and memory leaks to resolve before merging.

---

Bugs

1. Wrong variable in spline deallocation — equil/lar.f:~102 ⚠️ Critical

IF(.NOT. sq_in%allocated)THEN
    CALL spline_alloc(sq_in,ma,4)
ELSE
    CALL spline_dealloc(sq)    ! BUG: should be sq_in
    CALL spline_alloc(sq_in,ma,4)
ENDIF

sq (the global equilibrium spline) is being deallocated instead of sq_in. This would corrupt all subsequent equilibrium evaluations if this code path is taken. Fix this →

2. Duplicate WRITE statement — gpec/gpout.f:1930-1931

DO ising=1,msing
    WRITE(out_unit,'(1x,f6.3,17(es17.8e3))')       ! line 1930: blank write
    WRITE(out_unit,'(1x,f6.3,17(es17.8e3))')       ! line 1931: actual write
$        singtype(ising)%q,singtype(ising)%psifac,
$        spots(ising), ...

Line 1930 is a spurious duplicate with no arguments — it writes a blank/garbage line before the actual data line. The first WRITE should be removed. Fix this →

3. Orphaned WRITE statement — gpec/gpout.f:1918

WRITE(out_unit,'(1x,a13,a8,1x,a12,I2)')    ! format but no arguments
WRITE(out_unit,'(1x,a12,1x,I4)')"msing =",msing

After removing jac_out/tmag_out from the continuation lines, the WRITE header was left behind. It now writes a malformed blank line. This WRITE statement should be removed entirely. Fix this →

4. out_ahg2msc warning fires before namelist is read — equil/equil.f:49

SUBROUTINE equil_loadnamelists(...)
    IF (out_ahg2msc) THEN         ! checked BEFORE equil.in is read
        WRITE(*,*) "WARNING: ..."
    ENDIF
    CALL ascii_open(in_unit,"equil.in","OLD")
    READ(UNIT=in_unit,NML=equil_output)   ! out_ahg2msc read here

The out_ahg2msc guard fires on the default value (.FALSE.) before equil.in is ever opened. If a user sets out_ahg2msc = .TRUE. in equil.in, the warning will never appear. Move the check to after CALL ascii_close(in_unit). Fix this →

---

Memory Leaks

5. localcoup_out and localcoup_out_vecs not deallocated — gpec/gpout.f:1063-1065

Both are allocated in the IF (osing<msing) block (lines 813-817) but the deallocation only covers localcoup_out_bvecs and localcoup_out_vals. localcoup_out and localcoup_out_vecs (and matmo) are leaked on every call to gpout_singcoup. Fix this →

6. bwpmns_rmatch not deallocated — gpec/gpout.f:3783

Allocated alongside bwpmns but the DEALLOCATE(xmns,ymns,xnomns,bnomns,bwpmns) at line 3783 does not include bwpmns_rmatch. Fix this →

7. cellinfos arrays never deallocated — rdcon/dcon_mod.f / rdcon/gal.f

Six cellinfos arrays are allocated in gal_alloc but are not freed in gal_dealloc. They're also written to rdcon_netcdf.f but never cleaned up.

---

Design/WIP Concerns

8. use_res_spot behavior is confusing (gpec/gpec.f:562-574)

When gal_flag=.TRUE. and use_res_spot=.FALSE. (the default), the code automatically forces use_res_spot=.TRUE. with a warning. This makes use_res_spot=f in gpec.in effectively a no-op when gal_flag=t. Either document this as intentional behavior and make the default use_res_spot=t, or allow the user to opt out by not forcing the override.

9. res_psilim_diff declared but never used (gpec/gpec.f:63, gpec/gpglobal.F:37)

Appears in the gpec_control namelist and the global variable list but is not referenced anywhere in the code. Remove it or implement it.

10. Commented-out boundary guards (gpec/gpout.f)

lpsi=respsi-spots(ising)/(nn*ABS(singtype(ising)%q1))
! IF (lpsi < psilow) lpsi=psilow     ! <-- commented out

These guards were presumably commented out for a reason (debugging?), but leaving them out risks out-of-bounds spline evaluations near the plasma boundary. Add a note explaining why they're disabled or re-enable them.

11. Deprecated rdcon_read now just calls gpec_stop (gpec/rdcon.f:30-60)

The old implementation is left behind as commented-out code, and the function body now only calls gpec_stop. If it's truly deprecated, the function (and the surrounding commented-out code) can be removed entirely to reduce clutter.

12. Duplicate "subprogram 2" numbering (equil/equil.f:26 and :85)

Both equil_loadnamelists and equil_read are labeled as "subprogram 2".

---

Positive Changes

• Per-surface sing_spots: Moving from a global sing_spot to a per-surface array with resistive layer width scaling (S^{-1/3} dependence) is scientifically correct.
• interp_np_res=60: Configurable high-resolution inner-layer grid greatly improves visualization. The tot_grids formula update correctly accounts for the resolution increase.
• Single-surface chirikov fix: The new hdist=MIN(respsi, 1-respsi)/2.0 for msing==1 prevents a potential segfault/undefined behavior.
• q1_r fix: Was being set to 0.0_r8 when MRE data unavailable; now uses actual shear from spline.
• galsol%q from RMATCH: Previously reconstructed from spline evals, now read directly from globalsol.bin — more consistent with the resistive solution grid.
• bicube_alloc/bicube_dealloc safety: The allocated flag and associated guards improve robustness.
• wvac*psio**2.0_r8: Correct type suffix prevents potential mixed-precision issue.
• Error messages in gal_lsode_int: Printing x0 and x1 before stopping helps diagnose grid issues.
• RMATCH eigenvalue passed to GPEC: Proper handling of complex eigenvalue for resistive coupling.

---

Summary

The PR has good physics motivation and several important bug fixes. Before merging, the critical sq vs sq_in deallocation bug in equil/lar.f and the duplicate WRITE statements in gpec/gpout.f should be fixed. The memory leaks and the out_ahg2msc ordering bug should also be addressed. The res_psilim_diff dead variable and the use_res_spot override logic should be clarified.
| Branch: RDCON_singcoup

[logan-nc]

Independent Review & Bug Fixes

I performed an independent review of this PR (cross-checking the Claude review from Mar 11) and found several confirmed bugs plus two additional memory leaks not previously identified. All fixes have been pushed as individual commits.

Bug Fixes

equil/lar.f — CALL spline_dealloc(sq) → CALL spline_dealloc(sq_in) (commit a88910c)
In the sq_in re-allocation branch, the global equilibrium spline sq was being destroyed instead of the local input spline. This would corrupt global state on any re-entrant call to LAR.

equil/equil.f — Move out_ahg2msc check to after namelist read (commit 28f6df6)
The deprecation warning and vac_memory assignment were evaluated before equil.in was read, so the check always saw the default value (.FALSE.). Users setting out_ahg2msc=.TRUE. in their namelist would get no warning and vac_memory would always be .TRUE. regardless.

equil/direct.f — Remove spurious ALLOCATE/DEALLOCATE of y_out (commit 5de3b4a)
y_out was correctly declared as an automatic array sized by enstep. The ALLOCATE/DEALLOCATE added in this PR used the undefined symbol nstep (from lar.f, not in scope here) and failed to mark the array ALLOCATABLE, causing a compilation error. Reverted to the develop branch pattern.

gpec/gpout.f — localcoup_out replaced by localcoup_out_vecs (commit 93c4704)
localcoup_out was allocated but never populated with data, causing uninitialized memory to be written to the cl_id NetCDF variable. This appears to be a naming artifact: localcoup_out was an older name that was updated to localcoup_out_vecs (which is correctly filled), but the old variable was not fully removed. Fixed by removing localcoup_out entirely and reading localcoup_out_vecs directly. Also fixed associated memory leaks: matms, matmo, and localcoup_out_vecs were allocated but missing from their DEALLOCATE statements.

gpec/gpout.f — Remove orphaned WRITE statements in gpout_singfld ASCII output (commit c370b14)
Two spurious WRITE calls left after a prior refactor: a header write with format '(1x,a13,a8,1x,a12,I2)' but no data arguments, and a duplicate data write in the ising loop with no I/O list that emitted a blank line before each valid data line.

gpec/gpout.f — Add bwpmns_rmatch to DEALLOCATE in gpout_bparmns (commit 5601c3f)
bwpmns_rmatch was allocated alongside bwpmns but missing from the DEALLOCATE statement.

rdcon/gal.f — Deallocate cellinfos arrays in gal_dealloc (commit 45c2c86)
Six cellinfos component arrays (icell, iintvl, etypes, etypes_int, x1, x2) allocated in gal_alloc were never freed in gal_dealloc.

install/DEPENDENCIES.inc — Add msing.o as prerequisite of deltac.o (commit 545b00f)
deltac.f uses msing_bounds from msing_mod, but msing.o was not listed as a build dependency. Wildcard source discovery produces alphabetical order, so deltac.f was compiled before msing.f, failing with "Symbol msing_bounds not found in module".

Notes on Previous Review

Two items from the Mar 11 Claude review were not reproduced:

• res_psilim_diff declared but unused — it IS used at gpec/idcon.f:290
• use_res_spot silently overridden — the override at gpec/gpec.f:569 does print a warning message, so it is not silent (though it is a design concern)

[logan-nc]

@matt-pharr I am happy if you are. Please check my bug fix commits before merging. In particular, I restored the qlow and qhigh functionality under the assumption it was removed as part of development checks and not intended to stay. But if there is something broken about these in RDCON now, then please add printed warnings and/or errors if/when users try to use them.

[logan-nc]

@matt-pharr you have a green check ✅ from the github checks and an approved review from me. Ready to hit merge?!