[d-burg]

Summary

This PR adds tearing-mode growth rates to SLAYER. For each rational surface STRIDE finds, SLAYER builds the layer parameters from STRIDE's output and the kinetic profiles. It then scans the complex frequency plane adaptively for zeros of the matching condition between the inner-layer response and the outer Δ'. This is done surface by surface, or for several surfaces at once through the coupled determinant det(Δ' − Δ_crit − S^(1/3)Δ(Q)).

The latest round brings in the fixes found while porting this code to Julia (OpenFUSIONToolkit/GPEC), so the Fortran is ready for v1.6. Each fix is its own commit with its measured effect.

Physics fixes ported from the Julia code

Fix	Julia PR	Effect on the DIII-D-like example
Coupled determinant evaluated surface k at Q·τ_ref/τ_k; correct is Q·τ_k/τ_ref	#463	Every non-reference branch was at the wrong frequency (3/1 at Re Q = −0.51 instead of −1.06)
Δ' converted from STRIDE's ψ_N reference length to the r_s reference of the layer: K^(1/2+α_i)·Δ'_ij·K^(α_j−1/2), K = r_s dψ_N/dr, α = √(−D_I)	#403	2/1 diagonal × 0.80; γ(2/1) 107 → 86 /s
D_norm used √(τ/(1+τ)) (T_i, T_e transposed); now √ι_e with ι_e = Q_e/(Q_e − Q_i)	#434	D_norm(2/1) 3.24 → 2.75; γ 86 → 79.5 /s
Signed shear in τ_H, W_d and χ‖ gave NaN on reverse-shear surfaces	#431	None at positive shear; reverse shear now gives the same result as forward shear
ω* lacked the toroidal mode number n	#433	None at n = 1; Q_e, Q_i correct at n ≥ 2
Uncoupled matching dropped Im(Δ'_kk)	—	γ(2/1) +0.1% before today's equilibrium fixes on develop (#325, #329, #335), negligible after; the uncoupled result is now the diagonal limit of the coupled one

Unless noted, the numbers use the uncoupled scan with dc_type="none", so that they are comparable across fixes.

New in this round

• Root extraction. Before, gamma_sol was a placeholder zero and the root had to be found outside the code. Zeros and poles are now located from the AMR cells by winding number and refined by a secant iteration confined to the cell. The most unstable root is reported per surface as γ [1/s] and ω [rad/s]. All roots and poles go to slayer_output_n<n>.nc (Q_roots, gamma_roots, omega_roots, Q_poles, …).
• E×B Doppler shift in the coupled determinant (rotation_flag, default on), following #463. Each surface's layer is evaluated in its own E×B frame using the kinetic file's ω_E, so the coupled root is in the lab frame: the 2/1 moves from −10.5 to +31.0 krad/s, against Julia's 30.99 krad/s. scan_center sets the window center and defaults to the middle of the shifted branches.
• Default dc_type="rfitzp". On the example this damps the 2/1, 3/1 and 4/1, at −379, −645 and −2158 /s.

Robustness

• Compiler and determinant:
  • The branch did not compile under gfortran 15 (LSODE Jacobian argument); fixed.
  • The coupled determinant accepts any number of surfaces (LAPACK LU above rank 3) instead of stopping above 3.
• Mode numbers:
  • SLAYER takes n from the STRIDE run.
  • A local array had hidden the module-level nn used for the Alfvén time.
  • Namelist runs without est_gamma no longer stop with n = 0.
• Failure handling:
  • LSODE failures in the layer solve return NaN instead of an unconverged Δ.
  • A misspelled dc_type stops the run.
  • Coupling without read_eq stops with a message.
  • msing_max is clipped to the number of surfaces.
• Inputs, labels and output:
  • zeff and mu_i come from slayer.in instead of being hard-coded.
  • slayer.f reads every namelist variable develop reads (inlu, inpr_prof, QPscan2_flag, QDscan2_flag, Qratio), so develop-format and OMFIT-written inputs still parse. bal_flag is read from SLAYER_DIAGNOSE, as on develop, and is also accepted in SLAYER_OUTPUT.
  • inlu defaults to 1e8 as on develop. The branch had set it to 0, which zeroed the Lundquist number when params_flag=f.
  • Surfaces are labeled by m; q_rational holds the real q.
  • The estimate-only path no longer crashes.
• Precision: _r8 on every real literal in params.f; results move at about 1e-7.

Removed

• dc_type="toroidal". Its STRIDE geometric factor (Delta_geo) was not correct in the r_s reference: on the example it came out within 1% of rfitzp, where the validated derivation (#441) gives 22–60% above. The option now stops with a message. A port of #441 can follow once it merges there.
• The RMATCH Julia-benchmark patches (match_detgrid, per-surface Δ dump); this PR no longer touches rmatch/.

Verification

Run on the DIII-D-like example (docs/examples/DIIID_ideal_example, which now includes a slayer.in and a README section):

• Reference-surface invariance. With a real, diagonal Δ', the coupled roots equal each surface's uncoupled root rescaled by τ_1/τ_k to 7 digits. With rotation on, they equal it plus τ_1·n·Ω_E,k to 6 digits.
• n = 2. Q_e/τ_k at the q = 2 surface is exactly twice the n = 1 value.
• Other modes and inputs. Namelist runs at sval = ±1.5 give identical roots. The template input/slayer.in runs as shipped, as do the a10 example and the estimate-only mode.
• Build. gpec still builds and links against libslayer. DEPENDENCIES.inc matches gen_deps.py, and no code extends past column 72.

No solver tolerance changed. The legacy riccati used by GPEC's br_th keeps develop's tolerances.

Known differences and follow-ups

• Fortran and Julia growth rates differ on this example (2/1: 80 vs 192 /s without Δ_crit). The inputs differ:

  • The Fortran r_s is the flux-surface-averaged minor radius; Julia defaults to the midplane.
  • Z_eff is 2 here and 1 in the Julia deck.
  • The χ settings differ.

  The pieces that should agree do: the frequencies, and the rescale and Doppler identities.

• The legacy bal_flag torque balance in slayer.f can return NaN br_th (seen on the a10 example). This predates this branch: develop fixed the same issue only in gslayer.f.

• Δ' symmetrization (Symmetrize STRIDE Δ matrices #191). No longer needed, and Symmetrize STRIDE Δ matrices #191 is being closed. Off the diagonal, STRIDE's Δ' is already Hermitian up to a per-surface rescaling. Such a rescaling leaves det(Δ' − D) unchanged for any diagonal D, so symmetrizing cannot move a SLAYER root. The only part of Fitzpatrick's normalization that does change the roots is the per-surface reference length, which this PR applies (the K, α conversion above). Details are in Symmetrize STRIDE Δ matrices #191.

• Possible follow-ups: a regression/compare script for SLAYER, and a port of #441.

🤖 Generated with Claude Code

[logan-nc]

Thanks for making the PR @d-burg!
Don't be afraid to post figures (even block diagrams or initial figures showing what interpolate_slice_at_Q means) and questions directly here (you can go into the "Files Changed" tab and make in-line comments like "I did this to get blah. Is there a better way?", "need blah here - where do I get it?", etc.). Having things on record is helpful for years down the road when we are trying to remember what we did and why we made the decisions we made.

[d-burg]

Appearance of `slice` at a given Re(Q) = omega_ExB

[d-burg]

Sorting approach taken to match inner layer delta with outer layer delta in a way that should be robust to negative numbers. This is the sorting key which re-sorts the gammas array correspondingly, meaning the index of the match is true for both deltas and gamma

[logan-nc]

How's it going @d-burg? Have kin and sq interfaces set up yet?

[d-burg]

@logan-nc it's going! Having some trouble parsing exactly which dependencies and function/routine calls are needed where (with all the nested equilibrium objects). My main questions are:

1. Do I need to read in the vmat file to get msing? (This is how I've seen it done so far in GPEC)
2. Does the way I've written the kin and sq initializations look about right to you? and
3. Could you point me to a location where a netcdf of the form output by STRIDE is read by a different code? I've gotten a bit lost trying to figure out the deltaprime matrix read-in, and STRIDE deallocates the matrix I'd want before concluding its run

These questions all pertain to lines ~30-110 of the new layerinputs.f. Thanks!

[logan-nc]

You need read_kin from this module as well

[logan-nc]

I think the only equilibrium quantity the read_kin subroutine uses is chi1 (the normalization factor poloidal flux and psi_n). This is a public variable of the pentrc inputs module, so you should be able to set it more efficiently than having to do all this set_eq stuff with all the matrices. If that is hard for any reason (I think you might have to rename the variable like how dcon.F deals with overlapping names), then just add op_chi1 as an optional variable in the subroutine and use that if it is included in the inputs. See the progressbar in inputs.f for an example of how to use optional arguments.

[logan-nc]

I don't think you need the action matrices or to initialize_pentrc. You are not actually going to run pentrc. You literally just want to read that ascii table and make the kin spline.

[logan-nc]

You never assign anything to vmat_filename? Isn't it just blank then?

[logan-nc]

Searching the repo for vmat_filename I only see it in rmatch/msing.f. It seems to be a special interface needed to do the eigenfunction inner/outer layer matching from resistive DCON. Not something standard. Not something in stride or normal dcon. It also doesn't seem to contain sq, which is what you need, right? Seems like you might be barking up the wrong tree here?

[logan-nc]

@logan-nc it's going! Having some trouble parsing exactly which dependencies and function/routine calls are needed where (with all the nested equilibrium objects). My main questions are:

1. Do I need to read in the vmat file to get msing? (This is how I've seen it done so far in GPEC)
2. Does the way I've written the kin and sq initializations look about right to you? and
3. Could you point me to a location where a netcdf of the form output by STRIDE is read by a different code? I've gotten a bit lost trying to figure out the deltaprime matrix read-in, and STRIDE deallocates the matrix I'd want before concluding its run

These questions all pertain to lines ~30-110 of the new layerinputs.f. Thanks!

Thanks for checking in! Here are some quick replies

1. I don't think vmat is what you want it to be. It is a custom rdcon thing as far as I can tell. To read the equilibrium, you want to emulate how dcon does it:
   

   GPEC/dcon/dcon.F

   Lines 97 to 103 in 377c0d6

   	c-----------------------------------------------------------------------
   	c open output files, read, process, and diagnose equilibrium.
   	c-----------------------------------------------------------------------
   	CALL ascii_open(out_unit,"dcon.out","UNKNOWN")
   	CALL equil_read(out_unit)
   	CALL equil_out_global
   	CALL equil_out_qfind

   
   then a little lower down, it also does CALL sing_find, which is where msing gets set. It is all very annoying because each subroutine sets a mystery subset of "global" variables. It takes time to trace everything around. Sorry! Part of the fun I suppose.

2. See (1) for how to set sq by calling the appropriate equil subroutines. See inline comments regarding kin.

3. We don't have any netcdf reading in the code yet! I would love to see you lead us on this front. The top two hits of googling "read a netcdf file in fortran" look good to me. I am not an expert and just copied examples to do the writing part :)

[logan-nc]

Thanks Daniel, this is a big step forward for SLAYER. I checked the physics and ran the code against develop. The new growth-rate machinery holds up wherever I could derive it. A few small items are left before this can go into 1.6.0.

What I checked

• Physics fixes from the PR table: all are implemented correctly:
  • the Q·τ_k/τ_ref scaling in the coupled determinant, with the Doppler offset applied after the rescale;
  • the Δ' conversion, which I re-derived from the GGJ |x|^(1/2±α) asymptotics;
  • D_I and α, which are consistent with mercier.f and sing.f;
  • the shear chain rule;
  • D_norm = D_fac·√(Q_e/(Q_e−Q_i)), which reduces to √(5/3)ρ_s;
  • ω* with n, matching gpout.f and sign_conventions.rst;
  • Im(Δ'_kk) kept in the uncoupled matching.
• The E×B frame is internally consistent. It matches the Julia value of 30.99 krad/s.
• No root found: the code returns NaN with a warning, not a silent zero.
• Merge and existing outputs:
  • develop merges cleanly;
  • it builds with gfortran;
  • for users who don't run SLAYER, all text outputs of the four solovev regression chains, DIIID and a10 are identical to develop, and existing NetCDF variables are bit-identical.
• The DIIID example reproduces the README growth rates to within 0.4%. Develop's input/slayer.in and the a10 slayer.in still parse and run.

Pushed to this branch (3aabb99)

These were small enough to fix directly:

• Build fix: the slayer version.inc rule now uses git describe --tags --match 'v[0-9]*', like the other codes. On a full clone the tps-pin-* tags produced a version string past column 72, and slayer_netcdf.f failed to compile. CI's shallow checkout has no tags, so it didn't catch this.

• lnLamb initialized to 17 (sglobal.f). gpec_slayer never calls params(), so with slayer_threshold_flag in GPEC it ran with lnΛ=0 and got infinite S and Q.

• The estimate loop now sets tau and P_tor per surface, as the matching loop does. Before, every surface used the last surface's values. On DIIID with est_gamma_flag=t, the 2/1–5/1 estimates change by 12–25%.

• params() now stops with a message that names the bad value and the surface, instead of producing NaN or Inf, when:

  • n_e, T_e or T_i is not positive;
  • the shear is zero;
  • chi_perp is not positive where the critical Δ' needs it.

  Before, the chi_perp case ended in an unexplained "Wd iteration did not converge".

Still to fix

1. The old namelist defaults are zeroed (slayer.f:220-275).

   • Develop had nn=1, inpr=1, inds=6, n_e=1e19, t_e=t_i=1e3, inQ=20, inQ_e=2, inQ_i=-2.6, inc_beta=0.7, intau=1, Q0=4, delta_n_p=(1e-2,1e-2), jnum=500.
   • An old input that omits these now runs on zeros: inpr=0 gives LSODE atol=0, and inds=0 gives 0/0 in the riccati boundary condition.
   • Please restore the develop defaults for the variables that already existed.

2. gamma_fac is never assigned in namelist mode (slayer.f:269), so the estimate is always 0. Fix: gamma_fac = rs*REAL(delta_prime)/tau_r after params().

3. The br_th_flag path doesn't use the new per-surface inputs consistently:

   • it uses the namelist scalar sval for every surface (slayer.f:1250);
   • it uses Q0=Q_e where the torque balance needs the E×B Q (:1235);
   • it passes the Fitzpatrick D_norm into the Cole ds slot of riccati() (:1231-1233, and the stabscan branch at :1004).

   For 1.6, either store shear, ds and Q per surface, or stop with a clear message when br_th_flag is combined with read_eq.

4. dc_type='lar' looks dimensionally off (params.f, the lar case). It works out to rfitzp × 1/r_s². Please fix it or remove the option.

5. Two of the fixes reached only the read_eq path.

   • gslayer.f:172 still uses the signed sval in τ_H, so reverse-shear surfaces give NaN there.
   • gslayer.f:181-183 and params.f:173-174 still lack n in Q_e/Q_i.

   These should get the same ABS() and n changes.

6. The README should cover the coupled scan. With the shipped msing_max=3 and scan_width=1.5, the 2/1 falls outside the window when coupling_flag=t. The README should say to use msing_max=2 or a wider scan_width.

Not blocking

• The root finder can miss very weakly growing roots. The matching ODE is singular at g = −iQ_e, close to where a marginal tearing root sits. With the default refinement, a root with |γ̂| ≲ 0.002 can share a leaf cell with that pole, so the winding number is 0 and the root is missed. The result is NaN with a warning, not a wrong number, so this is fine for 1.6.

Once 1–6 are in I'm happy to approve for 1.6.0.

[logan-nc]

Approving for 1.6.0. Daniel is out of reach before the release, so I closed the six items myself in 2100323 (after merging develop in c3f9dda). Each follows a pattern already on this branch rather than a new design choice; please read them over when you are back, Daniel.

What changed

1. Unset namelist inputs — kept the branch's strict zero defaults (zero = unset, or "compute" for inQ, inQ_e, inQ_i), since restoring develop's silent values would collide with the "nonzero overrides params()" semantics in the new paths. Instead a run that needs an unset input stops and names it with its v1.5 default, e.g. a v1.5 input without inpr, inds, intau:

   ERROR: slayer.in leaves required inputs unset:
       inpr =  0.000E+00  (v1.5 default: 1.0)
       inds =  0.000E+00  (v1.5 default: 6.0)
       intau =  0.000E+00  (v1.5 default: 1.0)
     Since v1.6 SLAYER sets none of these silently. ...
   

   Dimensional inputs are checked when params() runs on the namelist, Cole inputs when the single-surface riccati() block runs. That block is now also skipped for br_th_flag, which computes its own deltas per surface like the growth-rate modes.
2. gamma_fac in namelist mode — rs*Re(delta_prime)/tau_r after params(), as build_inputs forms it. The template with delta_prime=(5,0) now estimates 102.7 Hz instead of 0.
3. br_th_flag per surface — slayer_inputs_type carries shear_arr, ds_arr, Qexb_arr; the Br loop uses each surface's shear, the Cole ds (not D_norm) and the ExB Q (not Q_e), like the other riccati() call sites. DIII-D with read_eq=t: finite thresholds on all five surfaces (1.0e-3, 5.6e-4, 7.0e-5, 5.0e-3, 2.1e-4). The ExB Q is Qconv*omega as in params() and gpec_slayer, without the n of the coupled Doppler shift; reconciling those for n ≥ 2 is a follow-up.
4. dc_type='lar' — removed, as 'toroidal' was. Algebraically lar = rfitzp/r_s: the same formula without its r_s normalization, so dimensionally inconsistent and, corrected, a duplicate. It stops with a message; the templates list "rfitzp" or "none".
5. |s| and n in the other paths — gslayer.f (GPEC's singthresh_flag) uses ABS(sval) in the Alfvén time and n in Q_e, Q_i (GPEC passes the n = 1 diamagnetic frequencies, the same expression build_inputs multiplies by n); params.f carries n into the namelist omega_e, omega_i. Identity at n = 1 with positive shear.
6. README — the coupled scan needs msing_max=2 or a wider scan_width; with msing_max=2 the coupled root is ω = 30.93 krad/s, γ = −379.6 /s, the value quoted in the PR body.

Checked

• Built with gfortran 15 on macOS; CI green on 2100323.
• DIII-D uncoupled run: screen output and slayer_output_n1.nc bit-identical before and after (2/1 root Q = −1.162293 − 0.04163i, γ = −379.42 /s, as in the README). The a10 slayer.in and develop's input/slayer.in give identical output before and after.
• DIII-D gpec with the shipped singthresh_flag=t runs to normal termination and writes singflx_crit_slayer and the torque-balance files.
• dc_type="lar" and under-specified v1.5 inputs stop with the messages above; the namelist estimate, and br_th_flag in both modes, run.

Merging with a merge commit so the per-fix history stays.

[logan-nc]

Follow-up to the approval above. An independent review of 2100323 found nothing merge-blocking and two small gaps in the unset-input report, fixed in b7fb8a4: qval is now reported (it feeds the namelist ω*), and the dimensional inputs are checked whenever params_flag=t, since that path calls params() on the namelist even with read_eq=t. Comments in gslayer.f/params.f now say plainly that GPEC hands gpec_slayer the bare diamagnetic frequencies and that GPEC's own nn multiplies them there (slayer.in is never read by GPEC). DIII-D, a10 and the develop template are unchanged by both commits.

Noted for after 1.6, not changed here:

• The ExB Q (Qconv*omega, from kin%f(5)) carries no n in params(), gpec_slayer and build_inputs, while Q_e, Q_i do and the coupled Doppler shift uses n*Omega_E. For n > 1 the torque balance 2*pr*(Q0 - Q) mixes the two. Identical at n = 1.
• The stabscan with Pperp_Ptor_flag=f (slayer.f, riccati(iing_coarse,Q_e,Q_i,P_perp,c_beta,D_norm,tau,pe,...)) still passes the Fitzpatrick D_norm in riccati()'s Cole ds slot, the pattern item 3 fixed in the Br loop.
• input_flag=t with params_flag=f still runs the single-surface baseline riccati() on the namelist Cole inputs, so inpr, inds, intau are required there although the per-surface results don't use them.
• delta.f documents riccati()'s 6th argument as "magnetic shear ds"; it is the Cole normalised ion sound radius.

[logan-nc]

Supersedes my follow-up comment above (which still lists items as open). An independent review of 2100323 found nothing merge-blocking; everything it raised is resolved on the branch, so 1.6.0 leaves no SLAYER items dangling.

b7fb8a4 — the unset-input report now covers qval (it feeds the namelist ω*) and applies whenever params_flag=t, since that path calls params() on the namelist even with read_eq=t. Comments in gslayer.f/params.f say plainly that GPEC hands gpec_slayer the bare frequencies and that GPEC's own nn multiplies them there (slayer.in is never read by GPEC).

Third commit —

• The normalised ExB Q carries n (Qconv*n*omega_E) in gpec_slayer and params(), as Q_e, Q_i and the coupled Doppler shift already did; before, the torque balance mixed an n-free Q0 with n-scaled Q_e for n ≥ 2. omega_sol is divided by n again. Identity at n = 1.
• The stability scan with Pperp_Ptor_flag=f gets the surface's Cole ds in riccati() instead of the Fitzpatrick D_norm, as the Br threshold does.
• The single-surface baseline riccati() is also skipped for input_flag, and the Cole inputs are required whenever a legacy scan will actually run.
• The standalone Br threshold (bal_flag, br_th_flag) ignores non-finite torque-balance entries and reports 0 with a warning instead of NaN, as gpec_slayer does. This closes the "standalone br_th can be NaN" known limitation in v1.6.0 release #350; the a10 example now prints br_th = 0 with the warning.
• MAXLOC counts from 1 while the scan arrays start at 0, so every reported Q_sol in slayer.f was one grid step high; fixed at all seven sites (a10's Q_sol moves from 0.05 to 0.0). gpec_slayer uses that point only as the centre of its refinement bracket, so it is unchanged and GPEC output is bit-identical.
• delta.f: riccati()'s inds is the normalised ion sound radius, not a shear.

Unchanged by all three commits: the DIII-D growth-rate run (screen and NetCDF) and develop's input/slayer.in. The two remaining SLAYER entries under Known limitations in #350 (delta_n_p hardcoded; the ~10× DIII-D threshold discrepancy) are a design choice and a physics investigation respectively, and stay documented there; any further work on them belongs in the Julia code.

[logan-nc]

@matt-pharr could you review and merge this one? Head is ba7a17f: develop merged in, my six review items closed in 2100323, and two follow-up commits (b7fb8a4, ba7a17f) from an independent re-review, all described in the comments above. CI is green on gfortran and ifort and my approval stands, but develop's ruleset wants a second approval because one of Daniel's older commits has an author email that isn't linked to a GitHub account. I'd rather have a second pair of eyes than an admin merge for the last Fortran release.

Where to look: slayer/slayer.f (the unset-input report after the namelist read, the Br-threshold loop, the MAXLOC subscript fix), slayer/gslayer.f (GPEC's path: ABS(sval), and GPEC's nn now multiplies the ExB and diamagnetic frequencies; identity at n = 1), slayer/params.f (same n, dc_type='lar' removed). Checked locally: the DIII-D growth-rate run and develop's input/slayer.in are bit-identical before and after; DIII-D gpec with singthresh_flag=t gives an identical gpec_singfld_n1.out; the a10 slayer.in goes from br_th = NaN to 0 with a warning and its Q_sol moves one grid step to the true maximum. Merge with a merge commit, please.

[matt-pharr]

Thank you @d-burg for this contribution! It is great to have your tearing workflow on hand in the public repo. I did not scrutinize your root finding logic but the changes to STRIDE are safe and everything seems to make relative sense and match project conventions. I made one tweak fixing the capitalization on the .out files, but that is it. Great work.