InnerLayer.SLAYER - FEATURE! - 🚨 Derive the toroidal critical-Δ geometric factor in the r_s reference

docs/development/references.md

@@ -108,3 +108,9 @@ The KineticForces module (formerly PENTRC) implements neoclassical toroidal visc
   - Published: Physical Review Letters **102**, 065002 (2009)
   - Link: https://doi.org/10.1103/PhysRevLett.102.065002
   - Describes: Trapped-particle nonambipolar transport theory underpinning the NTV calculation
+
+- **Connor et al. (2015)**: "The role of thermal conduction in tearing mode theory"
+  - Location: not in `docs/resources/`; preprint at https://arxiv.org/abs/1410.7240
+  - Published: Plasma Physics and Controlled Fusion **57**, 065001 (2015)
+  - Link: https://doi.org/10.1088/0741-3335/57/6/065001
+  - Describes: Toroidal critical Δ′ from finite parallel thermal conduction at a rational surface. Eq. 59 gives the geometric factor (`toroidal_dgeo`) and Eq. 20 the parallel-wavenumber gradient (`toroidal_kpar`) used by `dc_type = :toroidal` in `InnerLayer.SLAYER`; Eq. 61 is the large-aspect-ratio form of Eq. 59, and Eq. 65 gives the corresponding critical island width of Fitzpatrick (1995)

regression-harness/cases/diiid_slayer_n1.toml

@@ -109,6 +109,17 @@ label = "SLAYER Q_i"
 noise_threshold = 1e-12
 order = 21
 
+# Toroidal critical-Δ geometric factor, derived on every surface that carries a ResistGeometry
+# whatever dc_type is selected, so it is tracked here at the deck's dc_type = "none".
+[quantities.slayer_D_geo]
+h5path = "Tearing/PerSurface/D_geo"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER toroidal critical-Δ factor D_geo"
+noise_threshold = 1e-10
+order = 22
+class = "physics_converged"
+
 # Tearing eigenvalue (coupled mode → length 1). The headline deliverable:
 # a real, nonzero growth rate on a realistic equilibrium. Root-extraction
 # is sensitive to the AMR cell topology and ODE solver, so the Q/ω/γ

regression-harness/cases/diiid_slayer_n1_toroidal.toml

@@ -0,0 +1,70 @@
+# Regression case: DIII-D-like H-mode n=1 SLAYER with the toroidal critical-Δ offset.
+# Same deck as diiid_slayer_n1 with dc_type = "toroidal", so the geometric factor D_geo derived
+# from each surface's ResistGeometry feeds the χ_‖-matching offset D_c and, through it, the roots.
+# Each [quantities.*] block names an HDF5 path in the run output, how to extract it,
+# and the noise floor below which a difference is treated as zero.
+[case]
+name = "diiid_slayer_n1_toroidal"
+description = "DIII-D-like H-mode, n=1, SLAYER with the toroidal critical-Δ offset (dc_type = toroidal)"
+example_dir = "examples/DIIID-like_SLAYER_example"
+
+[overrides]
+"SLAYER.dc_type" = "toroidal"
+
+# Critical-Δ input and offset: the resistive interchange index and D_c. The geometric factor
+# D_geo does not depend on dc_type and is tracked once, in diiid_slayer_n1.
+[quantities.slayer_dr_val]
+h5path = "Tearing/PerSurface/D_R"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER resistive interchange D_R"
+noise_threshold = 1e-10
+order = 10
+class = "physics_converged"
+
+[quantities.slayer_D_c_offset]
+h5path = "Tearing/PerSurface/D_c_offset"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER critical-Δ offset D_c"
+noise_threshold = 1e-10
+order = 12
+class = "physics_converged"
+
+# Roots for the inner three rational surfaces, with the same floors as diiid_slayer_n1: the
+# outermost surfaces' contour search is unreliable, so they are not tracked.
+[quantities.slayer_Q]
+h5path = "Tearing/Roots/Q_root"
+type = "complex_vector"
+extract = "first_3_complex"
+label = "SLAYER Q_root [2/1,3/1,4/1]"
+noise_threshold = 1e-4
+order = 30
+class = "physics_converged"
+
+[quantities.slayer_gamma_Hz]
+h5path = "Tearing/Roots/gamma"
+type = "real_vector"
+extract = "first_3"
+label = "SLAYER γ_Hz [2/1,3/1,4/1]"
+noise_threshold = 1e-1
+order = 33
+class = "physics_converged"
+
+[quantities.slayer_no_root]
+h5path = "Tearing/Roots/no_root"
+type = "real_vector"
+extract = "first_3"
+label = "SLAYER no_root flags [2/1,3/1,4/1]"
+noise_threshold = 0
+order = 34
+class = "physics_converged"
+
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 999
+class = "diagnostic"

src/ForceFreeStates/Surfaces/ResistEval.jl

@@ -1,9 +1,9 @@
 # ResistEval.jl
 #
 # Per-singular-surface Glasser-Greene-Johnson geometric coefficients (E, F,
-# G, H, K, M) and the two flux-surface averages (⟨B²/|∇ψ|²⟩, ⟨B²⟩) that
-# downstream callers need to turn geometry into τ_A / τ_R with kinetic
-# profiles.
+# G, H, K, M) and the three flux-surface averages (⟨B²/|∇ψ|²⟩, ⟨B²⟩, ⟨|∇ψ|²⟩)
+# that downstream callers need to turn geometry into τ_A / τ_R with kinetic
+# profiles and into the toroidal critical-Δ geometric factor.
 #
 # Port of Fortran RDCON `resist_eval` (geometric part only).
 # Unlike the Fortran, this routine produces *only* the pure-equilibrium
@@ -19,6 +19,8 @@
 #   4: 1 / (B² · |∇ψ|²)
 #   5: B²
 #   6: |∇ψ|² / B²
+#   7: B          (see below)
+#   8: |∇ψ|²      (toroidal critical-Δ, ⟨|∇V|²⟩ = v1²⟨|∇ψ_N|²⟩)
 # All weighted by `jac / v1` (jacobian / dV/dψ) before integration.
 #
 # A seventh integrand, B, is added (beyond the Fortran set) so that ⟨B⟩ is
@@ -31,25 +33,27 @@
     ResistGeometry
 
 Per-singular-surface Glasser-Greene-Johnson geometric coefficients and
-supporting flux-surface averages.
-
-| field       | meaning                                              |
-|-------------|------------------------------------------------------|
-| `E`, `F`    | Glasser interchange parameters (enter `D_I = E+F+H-¼`) |
-| `G`         | Coupling coefficient (curvature × pressure gradient) |
-| `H`         | Pfirsch-Schlüter coefficient                         |
-| `K`         | Glasser parameter                                    |
-| `M`         | Mass factor                                          |
-| `avg_bsq_over_dpsisq` | ⟨B²/|∇ψ|²⟩ — needed for τ_R         |
-| `avg_bsq`   | ⟨B²⟩ — needed for τ_R                                |
-| `avg_B`     | ⟨B⟩ — needed for Lin-Liu-Miller f_t                  |
-| `B_max`, `B_min` | θ-extrema of B on the surface [T]               |
-| `f_trap`    | Lin-Liu & Miller 1995 trapped-particle fraction      |
-| `R_major`   | flux-surface-averaged major radius ⟨R⟩ [m]           |
-| `eps_local` | (R_max − R_min)/2 / R_major — local inverse aspect ratio |
-| `p_local`   | Plasma pressure at this surface [Pa]                 |
-| `p1_local`  | dp/dψ at this surface                                |
-| `v1_local`  | dV/dψ at this surface                                |
+supporting flux-surface averages. All averages ⟨·⟩ are normalized flux-surface
+averages (⟨1⟩ = 1), and `‖∇ψ_N‖` is the gradient of the normalized poloidal flux.
+
+| field                 | meaning                                                  |
+|:--------------------- |:-------------------------------------------------------- |
+| `E`, `F`              | Glasser interchange parameters (enter `D_I = E+F+H-¼`)   |
+| `G`                   | Coupling coefficient (curvature × pressure gradient)     |
+| `H`                   | Pfirsch-Schlüter coefficient                             |
+| `K`                   | Glasser parameter                                        |
+| `M`                   | Mass factor                                              |
+| `avg_bsq_over_dpsisq` | ⟨B²/‖∇ψ_N‖²⟩ — needed for τ_R                            |
+| `avg_bsq`             | ⟨B²⟩ — needed for τ_R                                    |
+| `avg_dpsisq`          | ⟨‖∇ψ_N‖²⟩ — needed for the toroidal critical-Δ           |
+| `avg_B`               | ⟨B⟩ — needed for Lin-Liu-Miller f_t                      |
+| `B_max`, `B_min`      | θ-extrema of B on the surface [T]                        |
+| `f_trap`              | Lin-Liu & Miller 1995 trapped-particle fraction          |
+| `R_major`             | flux-surface-averaged major radius ⟨R⟩ [m]               |
+| `eps_local`           | (R_max − R_min)/2 / R_major — local inverse aspect ratio |
+| `p_local`             | Plasma pressure at this surface [Pa]                     |
+| `p1_local`            | dp/dψ at this surface                                    |
+| `v1_local`            | dV/dψ at this surface                                    |
 
 `H` here is identical to the `H` reported by `mercier_scan!` and stored
 in `LocalStability/h` — the GGJ routine recomputes it for convenience.
@@ -69,6 +73,7 @@ struct ResistGeometry
     M::Float64
     avg_bsq_over_dpsisq::Float64
     avg_bsq::Float64
+    avg_dpsisq::Float64
     avg_B::Float64
     B_max::Float64
     B_min::Float64
@@ -83,10 +88,10 @@ end
 """
     resist_geometry(equil, psifac, q1; gamma=5/3) -> ResistGeometry
 
-Port of Fortran RDCON `resist_eval` restricted to the
-pure-equilibrium geometric coefficients. Integrates the 6 theta integrands
-at the given flux surface and combines them into E, F, G, H, K, M via the
-standard GGJ formulas.
+Glasser, Greene & Johnson 1975 (Phys. Fluids 18, 875) resistive-layer coefficients,
+restricted to the pure-equilibrium geometric part. Integrates the 8 theta integrands
+at the given flux surface and combines the 6 GGJ ones into E, F, G, H, K, M via
+the standard GGJ formulas.
 
 # Arguments
 
@@ -99,6 +104,7 @@ standard GGJ formulas.
   - `gamma`  — adiabatic index (default 5/3)
 
 !!! note "Contract"
+
     `psifac` must be a genuine interior rational surface (`0 < ψ < 1`) with
     nonzero `q1`, `p1 = dp/dψ`, and `p`. The GGJ combination divides by these
     and by `|∇ψ|²` (which → 0 at the axis), so calling on the magnetic axis,
@@ -107,47 +113,47 @@ standard GGJ formulas.
     rationals.
 """
 function resist_geometry(equil::Equilibrium.PlasmaEquilibrium,
-                          psifac::Real, q1::Real; gamma::Real=5/3)
+    psifac::Real, q1::Real; gamma::Real=5 / 3)
     profiles = equil.profiles
-    twopi    = 2π
-    chi1     = twopi * equil.psio
-    psi_f    = Float64(psifac)
+    twopi = 2π
+    chi1 = twopi * equil.psio
+    psi_f = Float64(psifac)
 
     # Surface-profile quantities (evaluate via the existing splines)
     twopif = profiles.F_spline(psi_f)
-    p      = profiles.P_spline(psi_f)
-    p1     = profiles.P_deriv(psi_f)
-    v1     = profiles.dVdpsi_spline(psi_f)
-    v2     = profiles.dVdpsi_deriv(psi_f)
-    q      = profiles.q_spline(psi_f)
-
-    # Build the 6 GGJ θ-integrands plus a 7th (B) for the neoclassical
-    # resistivity f_t calculation, and accumulate running extrema of
+    p = profiles.P_spline(psi_f)
+    p1 = profiles.P_deriv(psi_f)
+    v1 = profiles.dVdpsi_spline(psi_f)
+    v2 = profiles.dVdpsi_deriv(psi_f)
+    q = profiles.q_spline(psi_f)
+
+    # Build the 6 GGJ θ-integrands plus B (neoclassical resistivity f_t) and
+    # |∇ψ|² (toroidal critical-Δ factor), and accumulate running extrema of
     # (B, R) for Lin-Liu-Miller f_t and the local ε.
     ntheta = length(equil.rzphi_ys)
-    ff     = zeros(Float64, ntheta, 7)
-    B_max  = -Inf
-    B_min  =  Inf
-    R_max  = -Inf
-    R_min  =  Inf
+    ff = zeros(Float64, ntheta, 8)
+    B_max = -Inf
+    B_min = Inf
+    R_max = -Inf
+    R_min = Inf
     for itheta in 1:ntheta
         theta = equil.rzphi_ys[itheta]
-        f1  = equil.rzphi_rsquared((psi_f, theta))
-        f2  = equil.rzphi_offset((psi_f, theta))
+        f1 = equil.rzphi_rsquared((psi_f, theta))
+        f2 = equil.rzphi_offset((psi_f, theta))
         jac = equil.rzphi_jac((psi_f, theta))
         fy1 = FastInterpolations.deriv_view(equil.rzphi_rsquared, (0, 1))((psi_f, theta))
-        fy2 = FastInterpolations.deriv_view(equil.rzphi_offset,   (0, 1))((psi_f, theta))
-        fy3 = FastInterpolations.deriv_view(equil.rzphi_nu,       (0, 1))((psi_f, theta))
+        fy2 = FastInterpolations.deriv_view(equil.rzphi_offset, (0, 1))((psi_f, theta))
+        fy3 = FastInterpolations.deriv_view(equil.rzphi_nu, (0, 1))((psi_f, theta))
 
         rfac = sqrt(f1)
-        eta  = twopi * (theta + f2)
-        r    = equil.ro + rfac * cos(eta)
+        eta = twopi * (theta + f2)
+        r = equil.ro + rfac * cos(eta)
 
         v21 = fy1 / (2 * rfac * jac)
         v22 = (1 + fy2) * twopi * rfac / jac
         v23 = fy3 * r / jac
         v33 = twopi * r / jac
-        bsq    = chi1^2 * (v21^2 + v22^2 + (v23 + q*v33)^2)
+        bsq = chi1^2 * (v21^2 + v22^2 + (v23 + q * v33)^2)
         dpsisq = (twopi * r)^2 * (v21^2 + v22^2)
 
         B_here = sqrt(bsq)
@@ -163,6 +169,7 @@ function resist_geometry(equil::Equilibrium.PlasmaEquilibrium,
         ff[itheta, 5] = bsq
         ff[itheta, 6] = dpsisq / bsq
         ff[itheta, 7] = B_here
+        ff[itheta, 8] = dpsisq
         @views ff[itheta, :] .*= jac / v1
     end
     # Snap the repeated endpoint exactly equal to the start
@@ -182,7 +189,7 @@ function resist_geometry(equil::Equilibrium.PlasmaEquilibrium,
              (twopif * q1 * chi1 / avg[5] - v2)
     F_coef = (p1 * v1 / (q1 * chi1^2))^2 *
              (avg[1] * avg[3] + (twopif / chi1)^2 *
-              (avg[1] * avg[4] - avg[2]^2))
+                                (avg[1] * avg[4] - avg[2]^2))
     H_coef = twopif * p1 * v1 / (q1 * chi1^3) * (avg[2] - avg[1] / avg[5])
     M_coef = avg[1] *
              (avg[6] + (twopif / chi1)^2 * (avg[3] - 1.0 / avg[5]))
@@ -192,9 +199,9 @@ function resist_geometry(equil::Equilibrium.PlasmaEquilibrium,
 
     return ResistGeometry(
         E_coef, F_coef, G_coef, H_coef, K_coef, M_coef,
-        avg[1], avg[5],
+        avg[1], avg[5], avg[8],
         avg_B, B_max, B_min, f_trap, R_major, eps_local,
-        p, p1, v1,
+        p, p1, v1
     )
 end
 
@@ -206,8 +213,8 @@ Populate `sing.restype` for every `SingType` in `intr.sing` using
 filled.
 """
 function resist_eval_all!(intr::ForceFreeStatesInternal,
-                           equil::Equilibrium.PlasmaEquilibrium;
-                           gamma::Real=5/3)
+    equil::Equilibrium.PlasmaEquilibrium;
+    gamma::Real=5 / 3)
     for sing in intr.sing
         sing.restype === nothing || continue
         sing.restype = resist_geometry(equil, sing.psifac, sing.q1; gamma=gamma)

src/InnerLayer/InnerLayer.jl

@@ -25,7 +25,7 @@ import .GGJ: delta_convergence, solution_profile, asymptotic_profile, q4_surface
 
 import .SLAYER: SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear
 import .SLAYER: riccati_del_s, slayer_layer_thickness, LayerWidths
-import .SLAYER: surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs
+import .SLAYER: surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs, toroidal_dgeo, toroidal_kpar
 
 export InnerLayerModel, InnerLayerParameters, InnerLayerResponse, solve_inner, solve_inner_profile
 export GGJ, GGJModel, GGJParameters
@@ -37,6 +37,6 @@ export delta_convergence, solution_profile, asymptotic_profile, q4_surface_bench
 
 export SLAYER, SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear
 export riccati_del_s, slayer_layer_thickness, LayerWidths
-export surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs
+export surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs, toroidal_dgeo, toroidal_kpar
 
 end # module InnerLayer