diff --git a/docs/development/references.md b/docs/development/references.md index e19b92e32..43f3daf74 100644 --- a/docs/development/references.md +++ b/docs/development/references.md @@ -68,6 +68,19 @@ The PerturbedEquilibrium module implements GPEC-style perturbed equilibrium calc - Published: Physics of Plasmas **24**, 032505 (2017) - Describes: Self-consistent coupling with neoclassical effects +## InnerLayer Module (SLAYER) + +The SLAYER slab inner-layer solver and the resistive-layer width diagnostics are based on: + +- **Burgess et al. (2026)**: "Tearing Stability Prediction Combining Toroidal Calculations With a Two-Fluid Slab Layer" + - Location: `docs/resources/2026-Burgess-Tearing Stability Prediction Combining Toroidal Calculations With a Two-Fluid Slab Layer.pdf` + - Describes: The two-fluid slab layer, its four regimes, and the `S^(1/3)` outer-inner matching + +- **Fitzpatrick (2025)**: "Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation" + - Location: `docs/resources/2025-Fitzpatrick-Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation.pdf` + - Published: Nuclear Fusion (2025), doi 10.1088/1741-4326/ae4fdd (open access: arXiv:2511.07666) + - Describes: Diffusive-resistive layer width (Eq. 100) and the resistive-layer-overlap truncation `0 < Ψ < 1 − ε_c` (Sect. 5.9) + ## Resistive MHD Stability Analysis (Future Work) GPEC will eventually implement resistive MHD stability analysis based on: diff --git a/docs/resources/2025-Fitzpatrick-Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation.pdf b/docs/resources/2025-Fitzpatrick-Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation.pdf new file mode 100644 index 000000000..b12990f32 Binary files /dev/null and b/docs/resources/2025-Fitzpatrick-Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation.pdf differ diff --git a/docs/src/citations.md b/docs/src/citations.md index e78370553..dd8ef7546 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -131,6 +131,29 @@ Combines the toroidal outer-region calculation with the two-fluid slab layer mod --- +> R. Fitzpatrick, "Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation," +> *Nuclear Fusion* (2025). +> DOI: [10.1088/1741-4326/ae4fdd](https://doi.org/10.1088/1741-4326/ae4fdd) (open access: arXiv:2511.07666) + +Gives the diffusive-resistive layer width (Eq. 100) and, in Sect. 5.9, restricts the response calculation to the region inside the point where adjacent resonant layers first overlap. Both are implemented by the layer-overlap scan in `src/Tearing/LayerOverlap.jl`, which feeds the `psilim_from_layer_overlap` cap described in [Integration domain truncation](@ref). + +--- + +> R. Fitzpatrick, "Investigation of tearing mode stability near ideal stability boundaries via asymptotic matching techniques," +> *Physics of Plasmas* **32**, 062509 (2025). +> DOI: [10.1063/5.0273568](https://doi.org/10.1063/5.0273568) (arXiv:2503.24184) + +The critical island width `W_d` and its short/long-mean-free-path parallel conductivity blend, used for the critical-Δ offset in `src/InnerLayer/SLAYER/LayerParameters.jl`. + +--- + +> J. W. Connor, C. J. Ham, R. J. Hastie, and Y. Q. Liu, "The role of thermal conduction in tearing mode theory," +> *Plasma Physics and Controlled Fusion* **57**, 065001 (2015). arXiv:1410.7240 + +The toroidal (Eq. 59) and large-aspect-ratio (Eq. 61) critical-Δ offsets from finite parallel thermal conduction, building on H. Lütjens, J.-F. Luciani, and X. Garbet, *Physics of Plasmas* **8**, 4267 (2001). + +--- + ## ErrorFields Module > N. C. Logan et al., "Robustness of the tokamak error field correction tolerance scaling," diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index 858bb53a0..d5df873e0 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -12,8 +12,9 @@ extension: - **HDF5 (`.h5`/`.hdf5`)** — the GPEC kinetic schema (recommended). Datasets at the file root: required `psi` (normalized poloidal flux), `n_e`, `T_i`, `T_e`, `omega_E` (and `n_i`, defaulting to `n_e` if omitted); optional - `omega_tor`, `chi_e` (perpendicular heat diffusivity ``\chi_\perp``), and - `chi_phi` (toroidal momentum diffusivity ``\chi_\phi``). Each dataset + `omega_tor`, `chi_e` (perpendicular heat diffusivity ``\chi_\perp``), + `chi_phi` (toroidal momentum diffusivity ``\chi_\phi``), and `zeff` + (effective charge ``Z_\mathrm{eff}``, dimensionless). Each dataset carries a `units` attribute; the root carries `schema_version` and `provenance`. Densities are m⁻³, temperatures eV, frequencies rad/s, diffusivities m²/s. Additional datasets named `n_*` (e.g. `n_D`, `n_T`) @@ -24,9 +25,12 @@ extension: `psi_n n_i n_e T_i[eV] T_e[eV] omega_E`, retained for backward compatibility. Header rows are skipped. -The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_phi`, +The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_phi`/`zeff`, when present, are carried for the resistive-layer (SLAYER) analysis and ignored -here. +here. That analysis takes ``Z_\mathrm{eff}(\psi)`` from the `zeff` dataset when the file +has one, otherwise from `n_i/n_e` with the `[SLAYER]` impurity charge `zimp`, otherwise +from the scalar `[SLAYER] zeff`. Supply `zeff` whenever the impurity content is not a +single species fixed by `n_i/n_e`. ## Multi-ion runs diff --git a/docs/src/stability.md b/docs/src/stability.md index c61ca716e..4077ef68d 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -72,8 +72,8 @@ integrator = "forward" `riccati_eulerlagrange_integration` (the default) is our implementation of the STRIDE approach [Glasser 2018b], built on the dual Riccati reformulation [Glasser 2018a]. It -decomposes the radial domain into -independent chunks, integrates each chunk's fundamental-matrix (FM) propagator in parallel +decomposes the radial domain into independent chunks, integrates each chunk's +fundamental-matrix (FM) propagator in parallel using `Threads.@threads`, then multiplies the propagators in order and applies each singular-surface crossing serially. It is the only driver that produces the inter-surface ``\Delta'`` matrix. Because chunk endpoints are all it stores, `u_store` is sparse and stays @@ -253,6 +253,111 @@ The Galerkin integrator computes the same quantity in the same PEST-3 convention publishes it on the same path, so downstream consumers (SLAYER among them) never branch on which formalism ran. +## Integration domain truncation + +The Euler–Lagrange integration runs from the axis out to `psilim`, which is set from the +target safety factor `qlim`: + +- `qlim` starts at `min(qmax, qhigh)`, where `qmax` is q at the equilibrium's `psihigh`. +- With `set_psilim_via_dmlim = true` (the default, for diverted equilibria where q → ∞ at the + separatrix), `qlim` moves to `(last rational + dmlim) / n`, so the integration ends a fixed + fraction of a rational spacing past the last surface instead of close to one. The offset is + then stepped back by `1/n` until it no longer exceeds the bound it started from, so `dmlim` + never extends the domain past `qmax`, `qhigh`, or the layer-overlap cap below. A bound that + falls on a rational, or less than `dmlim / n` above one, therefore excludes it: with + `dmlim = 0.2`, both `qhigh = 4.0` and `qhigh = 4.1` give `qlim = 3.2`. +- With `set_psilim_via_dmlim = false` (limited or analytic equilibria with finite edge q), + `qhigh` and `psihigh` set the domain directly. + +### Resistive-layer overlap + +Near the separatrix the rational surfaces crowd together faster than their resistive layers +narrow, so eventually adjacent layers overlap. Past that point no surface keeps a +well-separated inner region and the matched-asymptotic treatment is not defined; Sect. 5.9 of +Fitzpatrick (2025), listed under [InnerLayer Module](@ref) in the citations, restricts the +response calculation to the region inside it. + +GPEC locates that point with a scan of the rational surfaces of each toroidal mode number in +the run; for multi-`n` runs the innermost overlap point, over all `n`, is the one used. Each surface gets a layer solve from the kinetic profiles, and its +diffusive-resistive width (Fitzpatrick 2025, Eq. 100) is converted to a normalized-flux width +so it can be compared with the spacing between neighbouring surfaces. The cut is placed at the +last surface before two adjacent layers first touch. + +Surfaces beyond the equilibrium grid are located on the logarithmic divergence of q at the +separatrix, taken from the equilibrium's own poloidal flux rather than fitted. Each X-point on +the last closed flux surface contributes one term, + +```math +q(\psi) = q_K + \sum_i A_i \ln\frac{\psi_{X_i} - \psi_K}{\psi_{X_i} - \psi}, +\qquad A_i = \frac{F}{2\pi R_{X_i}\sqrt{-\det H_i}}, +``` + +where ``\psi_{X_i}`` is the flux at the saddle, ``H_i`` the Hessian of the poloidal flux there, +``F = R B_\phi``, and ``q_K`` the traced q on the outermost grid surface ``\psi_K``. One term +describes a single null and two a double null, connected or not. The X-points are found by +following each minimum of the poloidal field along the outermost surface to the saddle it leads +to; on a single-null shape this usually also reaches the far secondary X-point, whose term is +small. The law is used only if it reproduces the traced shear over the outer knots, and an +equilibrium with no X-point nearby is not extrapolated. A limited plasma that has an X-point +nearby is not told apart from a diverted one: one rational is placed past its boundary and the run +warns that no overlap could be established there, leaving the domain unchanged. Divertors with a +second null beside the primary one or on the leg (snowflake, x-point target) are not modelled. + +The X-point flux need not coincide with the equilibrium's boundary: a reconstruction that places +its boundary slightly inside the separatrix has finite q at ``\psi = 1``, and rational surfaces +between that boundary and the X-point lie in the real plasma but outside the equilibrium GPEC +holds. The scan places only the first of them, to judge whether the outermost surface inside is +clean, and scores it only when the kinetic profiles extend that far. `psihigh` remains the upper +limit throughout: the overlap point can only narrow the domain, and when no overlap can be +established the domain keeps its requested limit. + +TokaMaker's `save_eqdsk` writes such files when `lcfs_pad` (``p``) is set with the default +`truncate_eq=True`: the equilibrium is cut at ``1 - p`` and renormalised, so the X-point sits at +``\psi = 1/(1 - p)``. Writing with `truncate_eq=False` keeps the whole plasma; EFIT-written files +place their boundary on the X-point flux. + +The scan takes its kinetic profiles from the `[SLAYER]` `profile_file`, through the same control +and loader the SLAYER analysis uses, so both see the same profiles, the same HDF5 group and the +same χ⊥(ψ) — the layer widths cannot disagree between them. A deck without a `[SLAYER]` +`profile_file` does not get the scan. Because the file is resolved from the run's TOML, the scan +only runs on the `gpec.toml` entry point: the programmatic `solve(prob, alg)` path has no inputs to +resolve it from, so the cap cannot be applied there and setting the flag warns instead. Its result is always written to +`ForceFreeStates/LayerOverlap/`, but it only constrains the domain when +`psilim_from_layer_overlap = true`. It then caps `qlim` before the `dmlim` step, so `dmlim` +still selects the final surface from inside the cap. A cap beyond `psihigh` has no effect. + +The flag is off by default, so every shipped deck is unchanged by it. With the shipped +DIII-D-like profiles the overlap point is ψ ≈ 0.9985: outside the ideal and Riccati decks' +`psihigh = 0.995`, where the cap could not bind even if enabled, but inside the SLAYER deck's +0.9995, where enabling it would exclude the q = 7 surface. `examples/DIIID-like_truncation_example` +raises `psihigh` to 0.9999 with the flag on: the cap excludes q = 7 and q = 8, and the resulting +domain (`psilim ≈ 0.99498`, last surface q = 6) reproduces the hand-chosen `psihigh = 0.995` of the +ideal deck. Whenever `psilim` lands past the overlap point, with the flag on or off, GPEC warns; +on the shipped decks that includes the SLAYER deck. + +**Which perpendicular diffusivity.** The two formulas that take `chi_perp` want different +quantities. The layer width goes through `P_perp = τ_R/τ_E` with `τ_E = r²/χ_E`, and χ_E is the +one-fluid perpendicular *energy* diffusivity, the gradient-weighted combination of the electron and +ion channels, `χ_E = (n_e χ_e ∇T_e + n_i χ_i ∇T_i)/(n_e ∇T_e + n_i ∇T_i)`. The critical-Δ island +width `W_d ∝ (χ_⊥/χ_∥)^{1/4}` describes electron-temperature flattening against an electron parallel +conductivity, so its χ_⊥ is the electron channel. GPEC currently feeds both from one profile, the +kinetic file's `chi_e`, which in the shipped DIII-D-like data is the electron power-balance +diffusivity: the right channel for the critical Δ, and a proxy for χ_E in the layer width. This is +the same single-χ convention as Fitzpatrick's TJ code, which also assumes `T_e = T_i`. +Separating the two is deferred: it needs an ion `chi_i` in the kinetic-profile schema, and the +power-balance ion diffusivity available upstream is not yet usable (it comes back negative across +most of the radius on the cases checked). Since both widths scale as `χ^{1/4}`, a factor-2 error in +χ moves them by about 19%, but that can be enough to move the overlap cut by a surface. + +| Dataset | Contents | +|---|---| +| `psi`, `m`, `n`, `r_s` | scanned surfaces and their minor radius | +| `delta_s_abs`, `width_delta_s` | Riccati layer thickness, in metres and as a Δψ width | +| `width_visco`, `width_dr` | visco-resistive and diffusive-resistive (Eq. 100) Δψ widths; `width_dr` sets the cut | +| `extrapolated` | 1 where the surface lies beyond the equilibrium grid | +| `psilim_overlap`, `first_overlap_index` | the cut, and the first surface whose layer overlaps its neighbour | +| `applied` | 1 when the cap was passed to the truncation (`dmlim` / `qhigh` may still cut deeper) | + ## Configuration reference All `ForceFreeStates` options are set in the `[ForceFreeStates]` section of `gpec.toml`. diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index a8a99f4c7..3546d08b1 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -72,7 +72,8 @@ scan_mode = "amr" # Q-plane scan strategy: "amr" (adaptive re coupling_mode = "uncoupled" # "uncoupled" (per-surface) or "coupled" (multi-surface determinant) dc_type = "none" # Critical-Δ offset selector: "none", "lar", "rfitzp", or "toroidal" mu_i = 2.0 # Ion mass in proton-mass units (2.0 = deuterium) -zeff = 1.0 # Effective charge +zeff = 1.0 # fallback only; Z_eff(ψ) comes from the kinetic file (its zeff profile, else n_i/n_e with zimp) +zimp = 6.0 # Impurity charge closing quasi-neutrality when Z_eff(ψ) is derived from n_i/n_e (6.0 = carbon) chi_perp = 1.0 # fallback only; the kinetic file supplies χ⊥(ψ) chi_tor = 1.0 # fallback only; the kinetic file supplies χ_φ(ψ) pole_threshold_adaptive = true # SLAYER |Δ| spans many orders; adapt to 10·median(|Δ|) diff --git a/examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5 b/examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5 index 2f084f080..dab58dc3f 100644 Binary files a/examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5 and b/examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5 differ diff --git a/examples/DIIID-like_truncation_example/gpec.toml b/examples/DIIID-like_truncation_example/gpec.toml new file mode 100644 index 000000000..72ec72a70 --- /dev/null +++ b/examples/DIIID-like_truncation_example/gpec.toml @@ -0,0 +1,50 @@ +# DIII-D-like H-mode equilibrium — n=1 truncation-only case for the resistive-layer overlap cap. +# Reuses the equilibrium and kinetic profiles of the sibling DIIID-like_ideal_example (referenced +# by relative path; not duplicated). psihigh is pushed well past the layer-overlap point (psi ~ 0.9985, +# where adjacent resistive layers first touch), so psilim_from_layer_overlap binds: the cap pulls qlim +# to 7 and dmlim steps back to 6.2, excluding q=7 and q=8. Ideal free-boundary stability only, with no local +# stability, forcing or perturbed equilibrium: the energies are kept so a regression case can pin +# what the capped domain does to the physics, not just to the bookkeeping. + +[Equilibrium] +eq_filename = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file +eq_type = "efit" # Type of the input 2D equilibrium file +jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) +grid_type = "auto" # Radial grid packing type ("auto" = two-pass measured-curvature grid) +psilow = 1e-4 # Lower limit of normalized poloidal flux +psihigh = 0.9999 # Upper limit of normalized poloidal flux (past the psi ~ 0.9985 layer-overlap point, so the cap binds) +mpsi = 0 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) +psi_accuracy = 0.001 # Target relative accuracy of splined profile derivatives for the auto grid +mtheta = 256 # Number of poloidal grid points +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) +etol = 1e-10 # Error tolerance for equilibrium solver +force_termination = false # Terminate after equilibrium setup (skip stability calculations) + +[ForceFreeStates] +local_stability_flag = false # Perform local stability analysis (Mercier and ballooning) across the ψ profile +vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes +force_termination = true # Terminate after force-free states (skip perturbed equilibrium calculations) + +psiedge = 1.0 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) +qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) + +nn_low = 1 # Smallest toroidal mode number to include +nn_high = 1 # Largest toroidal mode number to include +delta_mlow = 8 # Expands lower bound of Fourier harmonics +delta_mhigh = 8 # Expands upper bound of Fourier harmonics + +eulerlagrange_tolerance = 1e-10 # Relative tolerance for ODE integration of Euler-Lagrange equations +singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced + +# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) +integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) +set_psilim_via_dmlim = true # Truncate at (last_rational_q + dmlim)/n — TRUE for diverted equilibria (q → ∞ at separatrix) +dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) +psilim_from_layer_overlap = true # Cap psilim where adjacent resistive layers overlap (the scan is recorded either way) + +[SLAYER] +# The growth-rate analysis is not run here; the section is present because the layer-overlap scan +# takes its kinetic profiles, and its chi_perp(psi), from the SLAYER profile_file. +enabled = false # Run the SLAYER tearing-mode analysis +profile_file = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5" # Kinetic-profile file (relative to the run dir), required when SLAYER is enabled diff --git a/regression-harness/cases/diiid_truncation_n1.toml b/regression-harness/cases/diiid_truncation_n1.toml new file mode 100644 index 000000000..94688b3ba --- /dev/null +++ b/regression-harness/cases/diiid_truncation_n1.toml @@ -0,0 +1,120 @@ +# Regression case: DIII-D-like n=1 truncation-only run with the resistive-layer overlap cap enabled. +# psihigh = 0.9999 puts the overlap point (psi ~ 0.99851) inside the domain, so the cap binds: qlim +# is pulled from 8.35 to 7.0 and dmlim then steps back to 6.2, excluding the q=7 and q=8 surfaces. This is the +# only case where psilim_from_layer_overlap changes the domain; every shipped deck leaves it off. Ideal free-boundary stability only (no local stability, forcing or PE). Tracks the domain +# selection, the scan's layer widths, and the energies on the capped domain, so a change in what the +# truncation does to the physics cannot pass unnoticed. +[case] +name = "diiid_truncation_n1" +description = "DIII-D-like, n=1, psihigh=0.9999 with the resistive-layer overlap cap binding (truncation only)" +example_dir = "examples/DIIID-like_truncation_example" + +# Layer-overlap scan output. applied = 1 means the overlap bound was active going into sing_lim!; +# psilim_overlap is the last clean surface before adjacent Eq. (100) layers first touch. +[quantities.overlap_applied] +h5path = "ForceFreeStates/LayerOverlap/applied" +type = "int_scalar" +extract = "value" +label = "layer-overlap cap applied" +noise_threshold = 0 +order = 10 + +[quantities.psilim_overlap] +h5path = "ForceFreeStates/LayerOverlap/psilim_overlap" +type = "real_scalar" +extract = "value" +label = "layer-overlap psilim" +noise_threshold = 1e-8 +order = 11 + +[quantities.first_overlap_index] +h5path = "ForceFreeStates/LayerOverlap/first_overlap_index" +type = "int_scalar" +extract = "value" +label = "first overlapping surface index" +noise_threshold = 0 +order = 12 + +[quantities.width_dr] +h5path = "ForceFreeStates/LayerOverlap/width_dr" +type = "real_vector" +extract = "all_real" +label = "diffusive-resistive layer widths (Δψ)" +noise_threshold = 1e-8 +order = 13 + +# Final integration domain after the cap and the dmlim offset. +[quantities.qlim] +h5path = "Info/qlim" +type = "real_scalar" +extract = "value" +label = "qlim" +noise_threshold = 1e-12 +order = 20 + +[quantities.psilim] +h5path = "Info/psilim" +type = "real_scalar" +extract = "value" +label = "psilim" +noise_threshold = 1e-10 +order = 21 + +[quantities.mpert] +h5path = "Info/mpert" +type = "int_scalar" +extract = "value" +label = "mpert" +noise_threshold = 0 +order = 22 + +# Singular surfaces retained inside the capped domain (q=8 must be excluded). +[quantities.msing] +h5path = "SingularSurfaces/rational_count" +type = "int_scalar" +extract = "value" +label = "# singular surfaces" +noise_threshold = 0 +order = 30 + +[quantities.sing_q] +h5path = "SingularSurfaces/rational_q" +type = "real_vector" +extract = "all_real" +label = "singular q values" +noise_threshold = 1e-8 +order = 31 + +[quantities.sing_psi] +h5path = "SingularSurfaces/rational_psi" +type = "real_vector" +extract = "all_real" +label = "singular psi locations" +noise_threshold = 1e-8 +order = 32 + +# Energies on the capped domain — what the truncation does to the physics, not just the bookkeeping. +[quantities.ep_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_plasma_energies" +type = "complex_vector" +extract = "real_first" +label = "plasma energy Re(ep[1])" +noise_threshold = 1e-10 +order = 40 + +[quantities.et_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" +type = "complex_vector" +extract = "real_first" +label = "total energy Re(et[1])" +noise_threshold = 1e-10 +order = 41 + +# Runtime +[quantities.runtime] +h5path = "" +type = "runtime" +extract = "value" +label = "Runtime (s)" +noise_threshold = 0.0 +order = 999 diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl index 9bb0bd3ae..25fdcc46d 100644 --- a/src/Equilibrium/Equilibrium.jl +++ b/src/Equilibrium/Equilibrium.jl @@ -3,6 +3,7 @@ module Equilibrium # --- Module-level Dependencies --- using Printf, OrdinaryDiffEq, DiffEqCallbacks, LinearAlgebra, HDF5 +using Statistics: median using Roots using TOML using FastInterpolations @@ -14,6 +15,7 @@ import ..Utilities include("EquilibriumTypes.jl") include("SolveTolerances.jl") include("GridRefinement.jl") +include("SeparatrixQLaw.jl") include("FluxSurfaceMetrics.jl") include("CoordinateInvariant.jl") include("ReadEquilibrium.jl") diff --git a/src/Equilibrium/GridRefinement.jl b/src/Equilibrium/GridRefinement.jl index e58eb81dc..b103fc3fe 100644 --- a/src/Equilibrium/GridRefinement.jl +++ b/src/Equilibrium/GridRefinement.jl @@ -57,6 +57,95 @@ const CORE_MODEL_PSI_MAX = 0.03 const EDGE_MODEL_PSI_MIN = 0.9 # θ-lines subsample stride for the 2D geometry channels const THETA_STRIDE = 8 +# --- edge q-law gate ---------------------------------------------------------------------------- +# Minimum knots in the edge band before a fit is attempted. +const EDGE_FIT_MIN_KNOTS = 4 +# Weak absolute backstop against a fit that describes nothing at all. It is deliberately not a +# discriminator: no absolute r2 separates diverted from limited edges (real reconstructions fit the +# log law at 0.929-0.977, the limited a10 deck at 0.972), so the discrimination is the relative +# r2_log > r2_linear comparison, and the pathology guard is the shear consistency below. +const EDGE_FIT_MIN_R2 = 0.5 +# `A` in q = q_bar + A*(ln(1-psi) - u_bar) is the local logarithmic shear, so the fit must satisfy +# A = -(1-psi)*q'(psi) on the knots it was fitted to. Accept a factor of this much disagreement +# between the fitted slope and the measured one; a fit driven by edge noise, or one describing an +# edge that does not actually diverge, misses by far more. +const EDGE_FIT_SHEAR_RATIO = 2.0 + +# Least-squares slope and coefficient of determination for y = a + b*x. +function _linfit_r2(x::Vector{Float64}, y::Vector{Float64}) + x_bar = sum(x) / length(x) + y_bar = sum(y) / length(y) + sxx = sum((x .- x_bar) .^ 2) + sxx > 0 || return (NaN, NaN, NaN, NaN) + b = sum((x .- x_bar) .* (y .- y_bar)) / sxx + ss_res = sum((y .- (y_bar .+ b .* (x .- x_bar))) .^ 2) + ss_tot = sum((y .- y_bar) .^ 2) + r2 = ss_tot > 0 ? 1 - ss_res / ss_tot : NaN + return (b, r2, x_bar, y_bar) +end + +""" + edge_q_law(equil; psi_max, psi_min=EDGE_MODEL_PSI_MIN, min_knots=EDGE_FIT_MIN_KNOTS, + min_r2=EDGE_FIT_MIN_R2, shear_ratio=EDGE_FIT_SHEAR_RATIO) + -> nothing | (; A, q_bar, u_bar, n_knots, r2_log, r2_linear) + +Least-squares fit of the separatrix edge law `q = q̄ + A·(ln(1−ψ) − ū)` over the equilibrium's +outer knots. It gates the grid-refinement edge density floor, which needs only to know whether q +diverges at the edge; it does not locate anything. Placing surfaces beyond the grid is +[`separatrix_q_law`](@ref), which takes the divergence from the X-point(s) instead of assuming +it sits at ψ = 1. + +Returns `nothing` when the diverging model does **not** describe this equilibrium's edge, so a +plasma with finite edge q is never extrapolated as if q blew up: + + - fewer than `min_knots` knots in the band; + - `A ≥ 0`, i.e. q not rising toward ψ = 1; + - `r2_log < min_r2` — a weak backstop; the log law describes nothing at all; + - `r2_log ≤ r2_linear` — a plain linear-in-ψ fit explains the edge q at least as well, which is + what a **limited** plasma looks like. This comparison carries no scale and is what separates + the shipped limited decks (Solovev 0.760 vs 0.9996 linear; LAR 0.904 vs 0.998) from the + diverted ones (DIII-D 0.996 vs 0.865, 0.999 vs 0.594). + +This is a test of the **model**, not a topology classification: it asks whether q diverges +logarithmically here, not whether an x-point exists. Geometric x-point detection is +[`classify_topology`](@ref), which is a separate concern. +""" +function edge_q_law(equil::PlasmaEquilibrium; + psi_max::Real=Float64(equil.profiles.xs[end]), + psi_min::Real=EDGE_MODEL_PSI_MIN, + min_knots::Int=EDGE_FIT_MIN_KNOTS, + min_r2::Real=EDGE_FIT_MIN_R2, + shear_ratio::Real=EDGE_FIT_SHEAR_RATIO) + xs = collect(Float64, equil.profiles.xs) + band = findall(x -> x >= psi_min && x < psi_max, xs) + if length(band) < min_knots + n_tail = max(min_knots, length(xs) ÷ 10) + band = filter(i -> xs[i] < psi_max, collect(max(1, length(xs) - n_tail + 1):length(xs))) + end + length(band) >= min_knots || return nothing + + q = [Float64(equil.profiles.q_spline(xs[i])) for i in band] + u = [log(1.0 - xs[i]) for i in band] + all(isfinite, u) && all(isfinite, q) || return nothing + + A, r2_log, u_bar, q_bar = _linfit_r2(u, q) + _, r2_linear, _, _ = _linfit_r2([xs[i] for i in band], q) + (isfinite(A) && A < 0) || return nothing # q must rise toward the edge + (isfinite(r2_log) && r2_log >= min_r2) || return nothing + (isfinite(r2_linear) && r2_log > r2_linear) || return nothing + + # Shear consistency: the fitted slope is the local logarithmic shear, so it must match the + # profile's own -(1-psi)*q' over the same knots. + hint = Ref(1) + A_meas = median([-(1.0 - xs[i]) * Float64(equil.profiles.q_deriv(xs[i]; hint=hint)) for i in band]) + isfinite(A_meas) && A_meas < 0 || return nothing + ratio = A / A_meas + (ratio >= inv(shear_ratio) && ratio <= shear_ratio) || return nothing + + return (A=A, q_bar=q_bar, u_bar=u_bar, n_knots=length(band), r2_log=r2_log, r2_linear=r2_linear) +end + + # Rational-surface bracketing (Δ′ robustness). The ideal-MHD Δ′ asymptotic matching samples the # cubic equilibrium splines' 2nd/3rd derivatives across each rational ψ_s over the matching stencil # [ψ_s − dpsi, ψ_s + dpsi], dpsi = singfac_min/|n·q′|. A cubic 3rd derivative is piecewise constant @@ -248,10 +337,15 @@ function _knot_density(equil::PlasmaEquilibrium; tau::Float64, kin::Union{Nothin # nodal data of the smallest flux surfaces is dominated by integration and axis # extrapolation error, so measured curvature is not trusted below the core split. dlog = (4.0 * tau)^(1 / 3) + # The edge floor encodes the DIVERGING edge law q ≈ -A·ln(1-ψ), so it is applied only where + # that model actually describes the equilibrium. A limited plasma has finite edge q and must + # not be packed as if q blew up. The density itself stays A-independent (that is the point of + # the form: uniform relative q′ error regardless of A) -- the fit supplies validity, not slope. + edge_diverges = edge_q_law(equil) !== nothing @inbounds for i in 1:n if xs[i] <= CORE_MODEL_PSI_MAX rho_s[i] = 1.0 / (dlog * xs[i]) - elseif xs[i] >= EDGE_MODEL_PSI_MIN + elseif edge_diverges && xs[i] >= EDGE_MODEL_PSI_MIN rho_s[i] = max(rho_s[i], 1.0 / (dlog * (1.0 - xs[i]))) end rho_s[i] = max(rho_s[i], 1.0 / H_TARGET_MAX) diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index 934256554..3d4132d72 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -66,8 +66,9 @@ temperatures eV, frequencies rad/s, diffusivities m²/s. | `T_i`/`T_e` | ion / electron temperature | | `omega_E` | ExB rotation | | `omega_tor` | toroidal rotation (optional) | -| `chi_e` | perpendicular heat diffusivity χ⊥ | +| `chi_e` | perpendicular heat diffusivity χ⊥ (shipped data: electron power-balance χ_e; see the stability docs) | | `chi_phi` | toroidal momentum diffusivity χ_φ | +| `zeff` | effective charge Z_eff(ψ) supplied directly (optional); the resistive-layer analysis prefers it to its closure | | `species_densities` | named per-species densities (e.g. `"n_D"`, `"n_T"`) for explicit multi-ion input; `nothing` for ASCII / no extras | | `provenance` | short string recording the source file/format | """ @@ -81,16 +82,17 @@ struct KineticProfileData omega_tor::Union{Nothing,Vector{Float64}} chi_e::Union{Nothing,Vector{Float64}} chi_phi::Union{Nothing,Vector{Float64}} + zeff::Union{Nothing,Vector{Float64}} species_densities::Union{Nothing,Dict{String,Vector{Float64}}} provenance::String end function KineticProfileData(; psi, n_i=nothing, n_e=nothing, T_i=nothing, T_e=nothing, - omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_phi=nothing, + omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_phi=nothing, zeff=nothing, species_densities=nothing, provenance="") _v(x) = x === nothing ? nothing : Float64.(collect(x)) return KineticProfileData(Float64.(collect(psi)), _v(n_i), _v(n_e), _v(T_i), _v(T_e), - _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_phi), species_densities, String(provenance)) + _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_phi), _v(zeff), species_densities, String(provenance)) end const _KINETIC_H5_EXTS = (".h5", ".hdf5", ".he5") @@ -98,7 +100,7 @@ const _KINETIC_H5_EXTS = (".h5", ".hdf5", ".he5") # Dataset units written into / expected from the HDF5 kinetic schema. const _KINETIC_H5_UNITS = Dict("psi" => "normalized poloidal flux", "n_i" => "m^-3", "n_e" => "m^-3", "T_i" => "eV", "T_e" => "eV", "omega_E" => "rad/s", - "omega_tor" => "rad/s", "chi_e" => "m^2/s", "chi_phi" => "m^2/s") + "omega_tor" => "rad/s", "chi_e" => "m^2/s", "chi_phi" => "m^2/s", "zeff" => "dimensionless") """ read_kinetic_file(path; group="/") -> KineticProfileData @@ -137,7 +139,7 @@ function _read_kinetic_h5(path::AbstractString; group::AbstractString="/") haskey(ats, "provenance") && (prov = string(read(ats["provenance"]))) # Non-standard datasets named `n_*` (e.g. "n_D", "n_T") are named per-species density # profiles for explicit multi-ion input; other names are reserved for future schema fields. - standard = ("psi", "n_i", "n_e", "T_i", "T_e", "omega_E", "omega_tor", "chi_e", "chi_phi") + standard = ("psi", "n_i", "n_e", "T_i", "T_e", "omega_E", "omega_tor", "chi_e", "chi_phi", "zeff") extras = Dict{String,Vector{Float64}}() for k in keys(g) (k in standard || !startswith(k, "n_")) && continue @@ -147,7 +149,7 @@ function _read_kinetic_h5(path::AbstractString; group::AbstractString="/") return KineticProfileData(; psi=Float64.(vec(read(g["psi"]))), n_i=rd("n_i"), n_e=rd("n_e"), T_i=rd("T_i"), T_e=rd("T_e"), omega_E=rd("omega_E"), omega_tor=rd("omega_tor"), - chi_e=rd("chi_e"), chi_phi=rd("chi_phi"), + chi_e=rd("chi_e"), chi_phi=rd("chi_phi"), zeff=rd("zeff"), species_densities=(isempty(extras) ? nothing : extras), provenance=prov) end end @@ -179,6 +181,7 @@ function write_kinetic_h5(path::AbstractString, data::KineticProfileData; put("omega_tor", data.omega_tor) put("chi_e", data.chi_e) put("chi_phi", data.chi_phi) + put("zeff", data.zeff) # Named per-species densities (multi-ion input); dataset names must be `n_*` to round-trip. if data.species_densities !== nothing for (name, v) in data.species_densities diff --git a/src/Equilibrium/SeparatrixQLaw.jl b/src/Equilibrium/SeparatrixQLaw.jl new file mode 100644 index 000000000..48df51aeb --- /dev/null +++ b/src/Equilibrium/SeparatrixQLaw.jl @@ -0,0 +1,172 @@ +# SeparatrixQLaw.jl +# +# Safety factor between the outermost grid surface and the separatrix of a diverted equilibrium, +# taken from the X-point(s) of the equilibrium's own poloidal flux rather than fitted. +# +# Near a separatrix a field line stalls poloidally at each X-point it passes, so +# +# q(ψ) = c + Σ_i A_i ln[1 / (ψ_Xi − ψ)], A_i = F / (2π R_Xi sqrt(−det H_i)), +# +# with ψ_Xi the flux at saddle i of ψ(R, Z), H_i its Hessian there, and F = R B_φ. The constant is +# fixed by the traced q on the outermost grid surface, so nothing is fitted. One term describes a +# single null, two coincident terms a connected double null, and a second term at larger flux a +# disconnected one. +# +# The saddles are selected from the last closed flux surface itself: each minimum of the poloidal +# field along the outermost traced surface is followed to the critical point it leads to. On a +# single-null shape this usually also reaches the far secondary X-point; its term is small. A +# limited plasma with an X-point nearby is not told apart from a diverted one. Divertors with a +# second null beside the primary one (snowflake) or on the leg (x-point target) are not modelled; +# the shear check below is the backstop. + +# Accepted disagreement between the law's dq/dψ and the traced profile's over the outer knots. A +# law missing one of two equal X-points is low by a factor of two, so the bound sits well inside it. +const SEPARATRIX_SHEAR_RATIO = 1.5 +# Outer knots the shear check is taken over. +const SEPARATRIX_CHECK_KNOTS = 4 +# Closest approach to the separatrix flux a surface is placed at, in normalized flux. +const SEPARATRIX_PSI_GAP = 1e-9 + +# Gradient and Hessian of the 2D flux at (r, z). +function _flux_derivatives(psi_in, r::Float64, z::Float64) + pr = Float64(psi_in((r, z); deriv=DerivOp(1, 0))) + pz = Float64(psi_in((r, z); deriv=DerivOp(0, 1))) + prr = Float64(psi_in((r, z); deriv=DerivOp(2, 0))) + prz = Float64(psi_in((r, z); deriv=DerivOp(1, 1))) + pzz = Float64(psi_in((r, z); deriv=DerivOp(0, 2))) + return (; pr, pz, prr, prz, pzz, det=prr * pzz - prz^2) +end + +# Follow |∇ψ|² downhill from a seed with damped Newton steps on ∇ψ = 0. Returns the critical point, +# or `nothing` when the descent leaves the flux grid or stalls short of one. +function _descend_to_critical_point(psi_in, ingest::DirectIngest, r::Float64, z::Float64, grad_tol::Float64) + for _ in 1:100 + d = _flux_derivatives(psi_in, r, z) + g = d.pr^2 + d.pz^2 + sqrt(g) <= grad_tol && return (; r, z, det=d.det) + d.det == 0.0 && return nothing + dr = -(d.pzz * d.pr - d.prz * d.pz) / d.det + dz = -(d.prr * d.pz - d.prz * d.pr) / d.det + t = 1.0 + stepped = false + while t >= 1e-6 + rn, zn = r + t * dr, z + t * dz + if ingest.rmin <= rn <= ingest.rmax && ingest.zmin <= zn <= ingest.zmax + dn = _flux_derivatives(psi_in, rn, zn) + if dn.pr^2 + dn.pz^2 < g + r, z, stepped = rn, zn, true + break + end + end + t /= 2 + end + stepped || return nothing + end + return nothing +end + +""" + separatrix_q_law(equil; psi_max=equil.profiles.xs[end], shear_ratio=SEPARATRIX_SHEAR_RATIO) + -> nothing | (; psi_x, A, r_x, z_x, psi_anchor, q_anchor, shear_ratio) + +Safety factor beyond the equilibrium's ψ grid from the X-point(s) on its last closed flux surface: + + q(ψ) = q_anchor + Σ_i A_i ln[(psi_x_i − psi_anchor) / (psi_x_i − ψ)] + +`psi_x` are the normalized fluxes of the selected saddles (ascending, so `psi_x[1]` is where q +diverges), `A` their logarithmic coefficients `F / (2π R_X sqrt(−det H))`, `r_x`/`z_x` their +positions, and the anchor is the traced q at `psi_max`. `psi_x[1]` need not equal 1: a +reconstruction may place its boundary slightly inside the X-point flux. + +Returns `nothing` when the law does not apply or cannot be trusted: + + - the equilibrium is not a direct (2D flux) one — an inverse equilibrium carries its own q; + - no minimum of the poloidal field on the outermost surface leads to a saddle outside the grid + (a plasma with no X-point nearby; a limited plasma that has one still gets a law); + - q is not rising at the anchor; + - the law's dq/dψ and the traced profile's differ by more than `shear_ratio` over the outer + knots (returned as `shear_ratio`, law over traced). +""" +function separatrix_q_law(equil::PlasmaEquilibrium; + psi_max::Real=Float64(equil.profiles.xs[end]), + shear_ratio::Real=SEPARATRIX_SHEAR_RATIO) + ingest = equil.ingest + ingest isa DirectIngest || return nothing + psi_in = cubic_interp((ingest.psi_xs, ingest.psi_ys), ingest.psi_rz; extrap=ExtendExtrap()) + xs = equil.profiles.xs + psi_anchor = min(Float64(psi_max), Float64(xs[end])) + q_anchor = Float64(equil.profiles.q_spline(psi_anchor)) + q_slope = Float64(equil.profiles.q_deriv(psi_anchor)) + (isfinite(q_anchor) && q_anchor * q_slope > 0) || return nothing + + # Poloidal field along the outermost traced surface. + psi_surface = Float64(equil.rzphi_xs[end]) + thetas = equil.rzphi_ys[1:(end-1)] + hint2d = (Ref(1), Ref(1)) + rs, zs, bp = similar(thetas), similar(thetas), similar(thetas) + for (i, theta) in enumerate(thetas) + r_minor = sqrt(equil.rzphi_rsquared((psi_surface, theta); hint=hint2d)) + theta_cyl = 2π * (theta + equil.rzphi_offset((psi_surface, theta); hint=hint2d)) + rs[i] = equil.ro + r_minor * cos(theta_cyl) + zs[i] = equil.zo + r_minor * sin(theta_cyl) + d = _flux_derivatives(psi_in, rs[i], zs[i]) + bp[i] = hypot(d.pr, d.pz) / rs[i] + end + + # Each minimum of the poloidal field is where the surface passes a null; follow it there. + n = length(bp) + grad_tol = 1e-10 * maximum(bp .* rs) + merge_dist = 1e-6 * (ingest.rmax - ingest.rmin) + f_edge = abs(Float64(ingest.sq_fs[end, 1])) + saddles = NamedTuple{(:psi, :A, :r, :z),NTuple{4,Float64}}[] + for i in 1:n + (bp[i] < bp[mod1(i - 1, n)] && bp[i] <= bp[mod1(i + 1, n)]) || continue + cp = _descend_to_critical_point(psi_in, ingest, rs[i], zs[i], grad_tol) + (cp === nothing || cp.det >= 0) && continue + psi_x = 1.0 - Float64(psi_in((cp.r, cp.z))) / ingest.psio + psi_x > psi_anchor || continue + any(s -> hypot(s.r - cp.r, s.z - cp.z) < merge_dist, saddles) && continue + push!(saddles, (psi=psi_x, A=f_edge / (2π * cp.r * sqrt(-cp.det)), r=cp.r, z=cp.z)) + end + isempty(saddles) && return nothing + sort!(saddles; by=s -> s.psi) + + law = (psi_x=[s.psi for s in saddles], A=[s.A for s in saddles], r_x=[s.r for s in saddles], + z_x=[s.z for s in saddles], psi_anchor=psi_anchor, q_anchor=q_anchor) + + # The law must reproduce the traced profile's own shear where the two overlap. + knots = filter(i -> xs[i] <= psi_anchor, eachindex(xs)) + isempty(knots) && return nothing + band = knots[max(1, end - SEPARATRIX_CHECK_KNOTS + 1):end] + hint = Ref(1) + ratio = median([separatrix_q_dqdpsi(law, xs[i]) / Float64(equil.profiles.q_deriv(xs[i]; hint=hint)) for i in band]) + (isfinite(ratio) && inv(shear_ratio) <= ratio <= shear_ratio) || return nothing + return merge(law, (; shear_ratio=ratio)) +end + +""" +dq/dψ on the separatrix law. +""" +separatrix_q_dqdpsi(law, psi::Real) = sign(law.q_anchor) * sum(A / (px - psi) for (A, px) in zip(law.A, law.psi_x)) + +""" +q on the separatrix law. +""" +separatrix_q(law, psi::Real) = + law.q_anchor + sign(law.q_anchor) * sum(A * log((px - law.psi_anchor) / (px - psi)) for (A, px) in zip(law.A, law.psi_x)) + +""" +ψ at which the separatrix law reaches `q_target`. Returns `nothing` when the target is not beyond +the anchor, or lies closer to the separatrix flux than `psi_gap`. +""" +function separatrix_q_psi(law, q_target::Real; psi_gap::Real=SEPARATRIX_PSI_GAP) + s = sign(law.q_anchor) + s * (q_target - law.q_anchor) > 0 || return nothing + psi_sep = law.psi_x[1] + psi_gap < psi_sep - law.psi_anchor || return nothing + # Solve in the logarithmic distance to the separatrix, where q is close to linear. + residual(u) = s * (separatrix_q(law, psi_sep - exp(u)) - q_target) + u_near, u_far = log(psi_gap), log(psi_sep - law.psi_anchor) + residual(u_near) > 0 || return nothing + return psi_sep - exp(find_zero(residual, (u_near, u_far), Roots.Brent())) +end diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index 5544b0e60..4dde78771 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -136,6 +136,7 @@ gpec.toml. - `numunorms_init::Int` - Initial array size for solution normalization data - `singfac_min::Float64` - Fractional distance from rational q at which ideal jump condition is enforced - `set_psilim_via_dmlim::Bool` - Truncate the integration domain at `(last_rational_q + dmlim) / n` rather than at `qhigh` / `psihigh`. Fortran STRIDE found that truncating ~20 % above the outermost rational (`dmlim = 0.2`) avoids a numerical kink instability in δW that appears when the integration ends too close to or just below a rational surface. **For diverted equilibria where q → ∞ at the separatrix** (e.g. DIII-D geqdsks, the bulk of production use) this costs negligible physical domain because rationals get arbitrarily dense near the LCFS — `set_psilim_via_dmlim = true` is the safe and recommended default. **For limited circular / analytical equilibria with finite q at the edge** (Solovev, LAR scans), rationals are sparse and 20 % above the last rational chops off too much edge, so set `set_psilim_via_dmlim = false` and let `qhigh` / `psihigh` control the truncation. Multi-`n` runs are not supported by this truncation (the "outermost rational + dmlim / n" depends on which `n`); when `set_psilim_via_dmlim = true` with `nn_low != nn_high`, `sing_lim!` warns and falls back to `qhigh` / `psihigh`. Default `true`. + - `psilim_from_layer_overlap::Bool` - Cap the integration domain where adjacent rational surfaces' resistive layers overlap, after Fitzpatrick, Nucl. Fusion 2025 Sect. 5.9: past that point no surface retains a well-separated inner region. Applied as an upper bound on `qlim` at the top of `sing_lim!`, so `dmlim` / `qhigh` still select the final surface from inside it. The scan runs whenever the run's `[SLAYER]` section names a readable `profile_file` — the same profiles and χ⊥(ψ) the SLAYER analysis uses — and is recorded under `ForceFreeStates/LayerOverlap/` either way; this flag only decides whether it constrains the domain. It resolves that file from the run's TOML, so the programmatic `solve(prob, alg)` entry point cannot run it and the flag is inert there (with a warning). Off on every shipped deck; with the shipped DIII-D-like profiles the overlap point (ψ ≈ 0.9985) lies outside the 0.995 decks' domain but inside the SLAYER deck's 0.9995. See the integration-domain truncation section of the stability documentation. Default `false`. - `dmlim::Float64` - Distance beyond last rational surface (normalised ∈ [0,1) in units of 1/n). Only used when `set_psilim_via_dmlim` is true. Fortran STRIDE convention is 0.2 (truncate 20 % of one rational-surface spacing above the last surface), retained here. - `sing_order::Int` - Order of singular layer (Frobenius) expansion at rational surfaces. Default 6 (Fortran STRIDE convention for Δ' calculations; lower values trade accuracy for speed). - `qhigh::Float64` - Integration terminated at q limit determined by minimum of qhigh and qa from equil @@ -173,6 +174,7 @@ gpec.toml. numunorms_init::Int = 100 singfac_min::Float64 = 1e-4 # Matches Fortran STRIDE; required nonzero for the Riccati path. set_psilim_via_dmlim::Bool = true # Safe default for diverted equilibria (most production use); set false for limited/analytical (LAR, Solovev). Auto-skipped for multi-n. See docstring. + psilim_from_layer_overlap::Bool = false # Cap psilim where adjacent resistive layers overlap. Opt-in; the scan is recorded either way. See docstring. dmlim::Float64 = 0.2 sing_order::Int = 6 qhigh::Float64 = 1e3 diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl index 18c45def7..2e8382760 100644 --- a/src/ForceFreeStates/Surfaces/Finding.jl +++ b/src/ForceFreeStates/Surfaces/Finding.jl @@ -93,7 +93,8 @@ function sing_find!(intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEqui end """ - sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium) + sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium; + psilim_cap=nothing) Compute and set integration ψ, q, and q' limits by handling cases where user truncates before the last singular surface. Performs a similar function to `sing_lim` @@ -101,16 +102,22 @@ in the Fortran code. Main differences include renaming of sas_flag -> set_psilim removing dW edge storage variables since we now store all integration terms in memory, and simplification of the logic. -The target value `qlim` is first determined from user-specified control parameters -(`ctrl.qhigh` or `ctrl.dmlim`), subject to the constraint that it does not exceed -`equil.params.qmax`. If `set_psilim_via_dmlim` is true, `qlim` is adjusted to the largest -rational surface such that `nq + dmlim < qmax`. If `qlim < qmax`, a Newton iteration is -performed to find the corresponding `psilim` to integrate to. +The target `qlim` starts at `min(equil.params.qmax, ctrl.qhigh)`. When `psilim_cap` is given +and lies inside the domain, `qlim` is further capped at `q(psilim_cap)`: the resistive-layer +overlap point of Fitzpatrick, Nucl. Fusion 2025 Sect. 5.9, past which no surface retains a +well-separated inner region. The cap only ever narrows the domain; one beyond `psihigh` is inert. + +If `set_psilim_via_dmlim` is true, `qlim` is then moved to `(last rational + dmlim) / n` and +stepped back by `1/n` until it no longer exceeds the bound above, so `dmlim` never extends the +domain past `qmax`, `qhigh` or the cap. A bound that falls exactly on, or less than `dmlim / n` +above, a rational excludes that rational. If `qlim < qmax`, a Newton iteration finds the +corresponding `psilim` to integrate to. Note that the Newton iteration will be triggered if either `set_psilim_via_dmlim` is true or `ctrl.qhigh < equil.params.qmax`. Otherwise, the equilibrium edge values are used. """ -function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium) +function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium; + psilim_cap::Union{Nothing,Real}=nothing) profiles = equil.profiles @@ -119,6 +126,16 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, intr.q1lim = profiles.q_deriv(profiles.xs[end]; hint=Ref(profiles.npts_minus_1)) intr.psilim = equil.params.psihigh_resolved + # Resistive-layer overlap caps qlim (not psilim), so dmlim below still picks a surface inside it. + if psilim_cap !== nothing && psilim_cap < intr.psilim + q_cap = profiles.q_spline(Float64(psilim_cap)) + if q_cap < intr.qlim + @info "Resistive-layer overlap caps the domain: qlim $(@sprintf("%.3f", intr.qlim)) -> " * + "$(@sprintf("%.3f", q_cap)) (psi $(@sprintf("%.6f", intr.psilim)) -> $(@sprintf("%.6f", Float64(psilim_cap))))" + intr.qlim = q_cap + end + end + # Optionally override qlim based on dmlim (Fortran sas_flag=t equivalent). The cutoff reads # the *resolved* toroidal range on `intr`, so callers must assign intr.nlow / intr.nhigh # before calling; an unresolved range is an error rather than a silent change of truncation @@ -134,10 +151,12 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, @info "Setting psilim via dmlim: initial qlim = $(@sprintf("%.3f", intr.qlim)), dmlim = $(@sprintf("%.3f", ctrl.dmlim))" # Normalize dmlim ∈ [0,1) dmlim = mod(ctrl.dmlim, 1.0) + # qlim already holds min(qmax, qhigh, overlap cap); dmlim must never step past any of them. + qlim_bound = intr.qlim intr.qlim = (trunc(Int, intr.nlow * intr.qlim) + dmlim) / intr.nlow - # Reduce qlim if above qmax - while intr.qlim > equil.params.qmax + # Step back one rational at a time until inside the bound, keeping the dmlim offset + while intr.qlim > qlim_bound intr.qlim -= 1.0 / intr.nlow end end diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 23de76dad..6322d2c79 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -253,7 +253,9 @@ function main_from_inputs( ffs_start = time() locstab, ballooning_boundary = run_local_stability(ctrl, equil) - metric, mats = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles; species=kf_species) + metric, mats, layer_overlap = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles; + species=kf_species, + overlap_inputs=_overlap_slayer_inputs(inputs, intr.dir_path)) ffs_result = run_force_free_states(ctrl, equil, mats, intr, metric; runtimes=runtimes) if ctrl.write_outputs_to_HDF5 @@ -263,7 +265,8 @@ function main_from_inputs( inputs=inputs, forcing_modes=forcing_modes_snapshot, locstab=locstab, - ballooning_boundary=ballooning_boundary + ballooning_boundary=ballooning_boundary, + layer_overlap=layer_overlap ) written_h5 = ctrl.HDF5_filename @info "Results written to $(ctrl.HDF5_filename)" @@ -541,8 +544,70 @@ function run_local_stability(ctrl::ForceFreeStatesControl, equil::Equilibrium.Pl return locstab, ballooning_boundary end +# Kinetic profiles for the overlap scan come from the [SLAYER] section, through the same control +# and loader the SLAYER analysis itself uses, so the scan sees the same profiles, the same HDF5 +# group and the same chi(psi). Returns `nothing` when there is no usable [SLAYER] profile_file. +function _overlap_slayer_inputs(inputs, dir_path::AbstractString) + (inputs isa AbstractDict && haskey(inputs, "SLAYER")) || return nothing + ctrl = Runner.slayer_control_from_toml(inputs["SLAYER"]) + isempty(ctrl.profile_file) && return nothing + loaded = try + Runner._load_profiles(ctrl, dir_path) + catch err + @warn "Layer-overlap scan skipped: could not read the [SLAYER] profile_file." exception = err + return nothing + end + # chi(psi) and Z_eff(psi) from the file where it carries them, else the control's scalar + # fallbacks -- the same precedence run_slayer uses. Every other plasma parameter the layer + # widths depend on is the run's own [SLAYER] setting, so the scan and the SLAYER stage see + # one plasma. + return (; profiles=loaded.profiles, + chi_perp=loaded.chi_perp === nothing ? ctrl.chi_perp : loaded.chi_perp, + chi_tor=loaded.chi_tor === nothing ? ctrl.chi_tor : loaded.chi_tor, + zeff=loaded.zeff === nothing ? ctrl.zeff : loaded.zeff, + bt=ctrl.bt, mu_i=ctrl.mu_i, theta=ctrl.theta_sample, + resistivity_model=Runner._build_resistivity_model(ctrl.resistivity_model), + lnLambda_form=ctrl.lnLambda_form) +end + +# The binding scan across toroidal mode numbers: the one with the innermost overlap point, since +# the domain must stop before the first overlap of any n. Scans that find no overlap do not bind; +# when none overlap, the first scan is returned so its surfaces are still recorded. +function _binding_overlap(scans) + isempty(scans) && return nothing + overlapping = filter(sc -> sc.psihigh !== nothing, scans) + isempty(overlapping) && return first(scans) + return overlapping[argmin([sc.psihigh for sc in overlapping])] +end + +# Run the resistive-layer overlap scan for each toroidal mode number of the run and return the +# binding one, or `nothing` when there are no [SLAYER] profiles. Physics refusals (a surface that +# cannot be scored, a limited edge) are handled inside the scan; what qualifies its verdict is +# logged as a warning and the per-surface detail as info. +function _layer_overlap_scan(overlap_inputs, intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEquilibrium) + overlap_inputs === nothing && return nothing + # Eq. (100) is not covariant -- it is anchored to the toroidal-flux label of the paper's + # Eq. (30), so the scan is driven in that label regardless of the SLAYER default. + scans = [ + Tearing.resistive_layer_overlap(equil, overlap_inputs.profiles; n_tor=n, rs_method=:flux, + chi_perp=overlap_inputs.chi_perp, chi_tor=overlap_inputs.chi_tor, + bt=overlap_inputs.bt, mu_i=overlap_inputs.mu_i, zeff=overlap_inputs.zeff, theta=overlap_inputs.theta, + resistivity_model=overlap_inputs.resistivity_model, lnLambda_form=overlap_inputs.lnLambda_form) + for n in intr.nlow:intr.nhigh + ] + for (n, scan) in zip(intr.nlow:intr.nhigh, scans) + for msg in scan.caveats + @warn "Layer-overlap scan (n = $n): $msg" + end + for msg in scan.notes + @info "Layer-overlap scan (n = $n): $msg" + end + end + return _binding_overlap(scans) +end + """ - prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles) -> (metric, mats) + prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles; overlap_inputs=nothing) -> (metric, mats, overlap) Set up the force-free-states solve on `intr`: integration limits, the surviving singular surfaces and their GGJ coefficients, the poloidal mode range, and the metric plus @@ -554,10 +619,31 @@ function prepare_force_free_states!( equil::Equilibrium.PlasmaEquilibrium, kf_ctrl::KineticForces.KineticForcesControl, kinetic_profiles; - species=nothing + species=nothing, + overlap_inputs=nothing ) + # Resistive-layer overlap: locate where adjacent rational surfaces' layers run into each + # other, and use it as an upper bound on the integration domain. The scan runs whenever + # kinetic profiles are readable so `ForceFreeStates/LayerOverlap/` always records the point, + # but it only constrains the domain when the user opts in. + overlap = _layer_overlap_scan(overlap_inputs, intr, equil) + psilim_cap = (ctrl.psilim_from_layer_overlap && overlap !== nothing) ? overlap.psihigh : nothing + if ctrl.psilim_from_layer_overlap && overlap === nothing + @warn "psilim_from_layer_overlap = true but no layer-overlap scan was available; the domain is untouched." + end + # Determine psilim and qlim (where we will integrate to) - sing_lim!(intr, ctrl, equil) + sing_lim!(intr, ctrl, equil; psilim_cap=psilim_cap) + + # Fires whether or not the cap was applied: past the overlap point no surface retains a + # well-separated inner region, so a domain reaching beyond it is worth flagging even when + # the user did not opt in to the cap. + if overlap !== nothing && overlap.psihigh !== nothing && intr.psilim > overlap.psihigh + @warn "Integration domain extends past the resistive-layer overlap point: psilim = " * + "$(@sprintf("%.6f", intr.psilim)) > $(@sprintf("%.6f", overlap.psihigh)). Adjacent resistive " * + "layers overlap there, so the matched-asymptotic treatment is not defined. " * + "Set psilim_from_layer_overlap = true to cap the domain at that point." + end # Find all singular surfaces in the equilibrium sing_find!(intr, equil) @@ -643,7 +729,7 @@ function prepare_force_free_states!( end end - return metric, mats + return metric, mats, overlap end """ @@ -829,7 +915,7 @@ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrat end locstab, ballooning_boundary = run_local_stability(ctrl, equil) - metric, mats = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, nothing) + metric, mats, _ = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, nothing) result = run_force_free_states(ctrl, equil, mats, intr, metric) if ctrl.write_outputs_to_HDF5 @@ -1388,7 +1474,8 @@ function write_outputs_to_HDF5( inputs::Union{Nothing,Dict{String,Any}}=nothing, forcing_modes::Union{Nothing,Vector{ForcingTerms.ForcingMode}}=nothing, locstab::Union{FastInterpolations.CubicSeriesInterpolant,Nothing}=nothing, - ballooning_boundary=(psi=Float64[], alpha=Float64[], alpha_critical=Float64[]) + ballooning_boundary=(psi=Float64[], alpha=Float64[], alpha_critical=Float64[]), + layer_overlap=nothing ) ctrl = result.control @@ -1510,6 +1597,30 @@ function write_outputs_to_HDF5( out_h5["LocalStability/alpha"] = ballooning_boundary.alpha out_h5["LocalStability/alpha_critical"] = ballooning_boundary.alpha_critical + # Resistive-layer overlap scan. Written whenever the scan ran, whether or not it + # constrained the domain, so a run always shows where layer physics would have cut. + if layer_overlap !== nothing + lo = layer_overlap + out_h5["ForceFreeStates/LayerOverlap/m"] = lo.m + out_h5["ForceFreeStates/LayerOverlap/n"] = lo.n + out_h5["ForceFreeStates/LayerOverlap/psi"] = lo.psi + out_h5["ForceFreeStates/LayerOverlap/r_s"] = lo.rs + out_h5["ForceFreeStates/LayerOverlap/delta_s_abs"] = lo.delta_s_m + out_h5["ForceFreeStates/LayerOverlap/width_delta_s"] = lo.width_delta_s + out_h5["ForceFreeStates/LayerOverlap/width_visco"] = lo.width_visco + out_h5["ForceFreeStates/LayerOverlap/width_dr"] = lo.width_dr + out_h5["ForceFreeStates/LayerOverlap/extrapolated"] = Int.(lo.extrapolated) + out_h5["ForceFreeStates/LayerOverlap/psilim_overlap"] = + lo.psihigh === nothing ? NaN : lo.psihigh + out_h5["ForceFreeStates/LayerOverlap/first_overlap_index"] = + lo.first_overlap === nothing ? -1 : lo.first_overlap + # "applied" records that the bound was active going into sing_lim!, not that it set + # the final psilim -- dmlim/qhigh may truncate deeper inside it. + out_h5["ForceFreeStates/LayerOverlap/applied"] = + Int(ctrl.psilim_from_layer_overlap && lo.psihigh !== nothing && + lo.psihigh < equil.params.psihigh_resolved) + end + # Write integration data: the ψ trace and integrator diagnostics from the raw ODE state, # the ξ profiles from the solution. Either may be absent (Galerkin has no ODE state; # Riccati has no ξ solution), in which case the datasets are written empty. diff --git a/src/HDF5Schema.jl b/src/HDF5Schema.jl index 20801702e..ee63a3812 100644 --- a/src/HDF5Schema.jl +++ b/src/HDF5Schema.jl @@ -136,6 +136,33 @@ const MAIN_H5_ANNOTATIONS = [ (; long_name="Glasser-Greene-Johnson resistive interchange criterion D_R", dims=("psi",), attach=(1 => "LocalStability/psi",)), "LocalStability/ballooning_Delta_prime" => (; long_name="high-n ballooning Δ' (distinct from the tearing Δ')", dims=("psi",), attach=(1 => "LocalStability/psi",)), + # Resistive-layer overlap scan (Fitzpatrick, Nucl. Fusion 2025, Sect. 5.9). Its own surface + # axis: the list includes surfaces extrapolated beyond the integration domain, so it does not + # align with the SingularSurfaces/ axis. + "ForceFreeStates/LayerOverlap/psi" => + (; long_name="rational-surface location scanned for resistive-layer overlap", scale="surface_overlap"), + "ForceFreeStates/LayerOverlap/m" => + (; long_name="poloidal mode number of each scanned surface", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)), + "ForceFreeStates/LayerOverlap/n" => + (; long_name="toroidal mode number of each scanned surface", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)), + "ForceFreeStates/LayerOverlap/r_s" => + (; long_name="minor radius of each scanned surface in the Fitzpatrick flux label", units="m", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)), + "ForceFreeStates/LayerOverlap/delta_s_abs" => + (; long_name="Riccati resistive layer thickness |δ_s|", units="m", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)), + "ForceFreeStates/LayerOverlap/width_delta_s" => + (; long_name="|δ_s| expressed as a normalized-flux width", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)), + "ForceFreeStates/LayerOverlap/width_visco" => + (; long_name="visco-resistive comparison scale as a normalized-flux width", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)), + "ForceFreeStates/LayerOverlap/width_dr" => + (; long_name="diffusive-resistive layer width (Fitzpatrick 2025 Eq. 100) as a normalized-flux width; the channel the criterion uses", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)), + "ForceFreeStates/LayerOverlap/extrapolated" => + (; long_name="1 where the surface was located beyond the equilibrium grid via the separatrix edge q-law", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)), + "ForceFreeStates/LayerOverlap/psilim_overlap" => + (; long_name="domain limit implied by the first layer overlap, at the last surface with a well-separated inner region (NaN when no overlap was found in range)"), + "ForceFreeStates/LayerOverlap/first_overlap_index" => + (; long_name="index into this group's surface axis of the first surface overlapping its inner neighbour (-1 when none do)"), + "ForceFreeStates/LayerOverlap/applied" => + (; long_name="1 when the overlap bound was handed to sing_lim! as an active cap on qlim (dmlim/qhigh may still truncate deeper); 0 when recorded only"), "LocalStability/ballooning_psi" => (; long_name="normalized poloidal flux ψ_N of the ballooning α boundary scan", scale="psi_ballooning"), "LocalStability/alpha" => diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 252e81ae9..87488e18e 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -83,13 +83,15 @@ derivatives are analytic — no label carries a finite-difference stencil. implementation divides it by 2π to form `g`. Defined from flux alone, it carries no circular-cross-section assumption, and its derivative `dr/dψ ∝ q` grows toward a separatrix where the - geometric labels' `da/dψ` collapses. + geometric labels' `da/dψ` collapses. Beyond the ψ grid the q spline is a + cubic extrapolation that does not follow that growth, so `q_beyond`, a + callable of ψ, supplies q there when the caller knows it. On shaped equilibria the labels agree at low-q surfaces and diverge strongly near the edge, where the slab-layer matching is label-ambiguous regardless of choice. The label is selected programmatically; it is not exposed via TOML. """ -function radial_label(equil; rs_method::Symbol=:midplane, theta::Real=0.0) +function radial_label(equil; rs_method::Symbol=:midplane, theta::Real=0.0, q_beyond=nothing) theta_f = Float64(theta) _flux_r, _flux_dr = if rs_method === :flux @@ -103,8 +105,9 @@ function radial_label(equil; rs_method::Symbol=:midplane, theta::Real=0.0) Phi = collect(Float64, cumulative_integrate(cubic_interp(xs_f, qg))) r_knots = sqrt.(max.(2 .* psiof .* Phi ./ b0f, 0.0)) rspl = cubic_interp(xs_f, r_knots; extrap=ExtendExtrap()) + _q_at(ψ) = (q_beyond !== nothing && ψ > xs_f[end]) ? Float64(q_beyond(ψ)) : Float64(equil.profiles.q_spline(ψ)) (ψ -> Float64(rspl(Float64(ψ))), - ψ -> psiof * Float64(equil.profiles.q_spline(Float64(ψ))) / + ψ -> psiof * _q_at(Float64(ψ)) / (b0f * _g_at(Float64(ψ)) * max(Float64(rspl(Float64(ψ))), eps()))) else (nothing, nothing) @@ -180,7 +183,8 @@ profiles, without an intermediate file round-trip. physical field, so passing it as a scalar is almost always wrong. - `mu_i` -- ion mass in proton-mass units (default `2.0` for D). - - `zeff` -- effective charge (default `1.0`). + - `zeff` -- effective charge. Scalar or a callable of `psi`, evaluated per + surface (default `1.0`). - `chi_perp` -- perpendicular heat diffusivity [m²/s]. Scalar or a callable of `psi` (default `1.0`). - `chi_tor` -- toroidal heat diffusivity [m²/s]. Scalar or a callable @@ -190,7 +194,7 @@ profiles, without an intermediate file round-trip. (`:lar`, `:rfitzp`, `:toroidal`). When `nothing` (default), Julia derives it per-surface from the equilibrium as `dr_val_k = D_R(ψ_k) = E_k + F_k + H_k²`, - consistent with Connor-Hastie-Helander 2015 (PPCF 57 065001) Eq. 59 + consistent with Connor, Ham, Hastie & Liu 2015 (PPCF 57 065001) Eq. 59 which uses `(−D_R)` in the χ_‖-matching critical-Δ. Pass a scalar / vector / callable to override. @@ -229,8 +233,9 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; bt=nothing, R0=nothing, rs_method::Symbol=:midplane, + q_beyond=nothing, mu_i::Real=2.0, - zeff::Real=1.0, + zeff=1.0, z_i::Real=1.0, chi_perp=1.0, chi_tor=1.0, @@ -255,7 +260,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; Float64(bt(ψ)) end - _rs_at, _da_dpsi_at = radial_label(equil; rs_method=rs_method, theta=theta) + _rs_at, _da_dpsi_at = radial_label(equil; rs_method=rs_method, theta=theta, q_beyond=q_beyond) # Per-surface ω_*e, ω_*i (diamagnetic frequencies) from spline derivatives. # In flux coordinates ω_* = n·(dp/dψ)/(e·n_e), so the surface's toroidal mode @@ -308,12 +313,12 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; else lnL = coulomb_log_e(prof.n_e, prof.T_e; form=lnLambda_form) nu_star_e(prof.n_e, prof.T_e, rg.R_major, rg.eps_local, - q, zeff; lnLamb=lnL) + q, _eval(zeff, psi); lnLamb=lnL) end # dr_val: per-surface resistive interchange index D_R = E + F + H² # (Glasser-Greene-Johnson 1975). Used by `_solve_dc_tmp` to compute - # the χ_‖-matching critical-Δ via Connor-Hastie-Helander 2015 Eq. 59, + # the χ_‖-matching critical-Δ via Connor et al. 2015 Eq. 59, # which has `(−D_R)` as a multiplier. NOT the Mercier index # D_I = E + F + H − 1/4 (see this function's docstring); we use the # physically correct D_R here. @@ -332,8 +337,8 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; _eval(dr_val, psi) end - # dgeo_val: only used by dc_type=:toroidal (the Connor-Hastie- - # Helander 2015 formula). Auto-derivation requires ⟨|∇ψ|²⟩ FSA + # dgeo_val: only used by dc_type=:toroidal (the Connor et al. + # 2015 formula). Auto-derivation requires ⟨|∇ψ|²⟩ FSA # which the current `ResistGeometry` doesn't expose; for now we # require an explicit value if the toroidal dc_type is selected. dgeo_val_k = if dgeo_val === nothing @@ -377,7 +382,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; n_e=prof.n_e, t_e=prof.T_e, t_i=prof.T_i, omega=prof.omega, omega_e=ω_e_use, omega_i=ω_i_use, qval=q, sval_r=sval_r, bt=_bt_at(psi), - rs=rs, R0=R0_use, mu_i=mu_i, zeff=zeff, + rs=rs, R0=R0_use, mu_i=mu_i, zeff=_eval(zeff, psi), chi_perp=_eval(chi_perp, psi), chi_tor=_eval(chi_tor, psi), m=m_res, n=n_res, diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl index ae875dea4..52ac0e173 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/src/InnerLayer/SLAYER/LayerParameters.jl @@ -41,7 +41,7 @@ de-normalization. The parametrization uses `P_perp`, `P_tor`, and | `bt` | Toroidal field [T] | | `sval_r` | r-based magnetic shear r_s · (dq/dr) / q (Fitzpatrick convention) | | `dr_val` | Resistive interchange D_R = E + F + H² (critical-Δ input; auto-derived from GGJ coefficients unless overridden) | -| `dgeo_val` | Connor-Hastie-Helander 2015 Eq. 59 geometric factor (0 unless supplied) | +| `dgeo_val` | Connor et al. (2015) Eq. 59 geometric factor (0 unless supplied) | | `eta` | Parallel resistivity entering τ_R = μ₀r_s²/η [Ω·m] | | `d_beta` | Beta-weighted ion length scale c_β · d_i [m] | | `dc_tmp` | Critical-Δ offset from chi_parallel matching | @@ -129,8 +129,10 @@ function r_based_shear(rs::Real, q::Real, dq_dpsi::Real, da_dpsi::Real) return rs * dq_dpsi / (q * da_dpsi) end -# Internal: solve the Wd self-consistency loop for the chi_parallel-based -# critical Δ (Connor-Hastie-Helander 2015). Returns dc_tmp as a Float64. +# Internal: solve the Wd self-consistency loop for the chi_parallel-based critical Δ. W_d and the +# smfp/lmfp chi_parallel blend follow Fitzpatrick, Phys. Plasmas 32, 062509 (2025); the toroidal and +# large-aspect-ratio offsets are Connor, Ham, Hastie & Liu, PPCF 57, 065001 (2015) Eqs. (59) and (61), +# after Lütjens, Luciani & Garbet, Phys. Plasmas 8, 4267 (2001). Returns dc_tmp as a Float64. function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real, chi_perp::Real, t_e::Real, zeff::Real, tau_ee::Real, rs::Real, R0::Real, sval_r::Real, n_tor::Integer, @@ -210,7 +212,9 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl - `R0` -- major radius [m] - `mu_i` -- ion mass in proton-mass units (e.g. 2.0 for D) - `zeff` -- effective charge - - `chi_perp`, `chi_tor` -- perpendicular / toroidal heat diffusivity [m²/s] + - `chi_perp`, `chi_tor` -- perpendicular / toroidal heat diffusivity [m²/s]. `chi_perp` feeds both + `P_perp` (which wants the one-fluid energy diffusivity χ_E) and the critical-Δ `W_d` (which wants + the electron channel); one profile serves both for now, see the stability docs - `m`, `n` -- poloidal / toroidal mode numbers at the surface - `dr_val`, `dgeo_val` -- inputs for the critical-Δ formula - `dc_type` -- one of `:none`, `:lar`, `:rfitzp`, `:toroidal` diff --git a/src/InnerLayer/SLAYER/LayerThickness.jl b/src/InnerLayer/SLAYER/LayerThickness.jl index 0b6db9761..a07d4ea48 100644 --- a/src/InnerLayer/SLAYER/LayerThickness.jl +++ b/src/InnerLayer/SLAYER/LayerThickness.jl @@ -137,6 +137,26 @@ is built from, retained as a drift-scale reference. - `delta_s` -- complex layer thickness `δ_s = dels_db · d_β` [m] - `delta_s_m` -- `|δ_s|`, the resistive layer thickness in meters (primary) - `d_beta` -- β-weighted ion scale `c_β·d_i` in meters (drift reference) + - `delta_norm` -- layer normalization length `r_s · S^(-1/3)` in meters, the + length by which Δ' is made dimensionless (`Δ̂' = Δ'·delta_norm`). Shared by all + four regimes of Burgess et al. (2026) Eqs. (10)-(13), not specific to any one + of them, and **not** the classic non-rotating Furth-Killeen-Rosenbluth (1963) + constant-ψ tearing width, which scales as `S^(-2/5)` and is not computed here. + Exactly the reciprocal of the `delta_n` Δ-normalization already carried on + `SLAYERParameters`, restated as a length for comparison with `delta_s_m` + - `delta_dr` -- diffusive-resistive layer thickness in meters, Fitzpatrick (2025) + Eq. (100). Strong magnetic shear near the separatrix forces every resonant layer in + that region into the diffusive-resistive (`nu = 1/4`) regime, so this is the width the + edge overlap criterion compares against surface spacing. In the paper's normalized + radius, `delta_hat = tau_A^(1/2) / (tau_R^(1/4) tau_E^(1/4) d_beta_hat^(1/2) (n|s|)^(1/2))`; + with `lu = tau_R/tau_A` and `P_perp = tau_R/tau_E` that is + `lu^(-1/2) P_perp^(1/4) / (d_beta_hat^(1/2) (n|s|)^(1/2))`, scaled by `rs` for meters. + Distinct from `delta_visco`, which is the 2/3 power of the same timescale grouping and + carries neither the `d_beta` nor the shear factor + - `delta_visco` -- viscous-resistive scale `delta_norm · P_perp^(1/6)` in meters, + from the `P^(1/6)` broadening of the VR-regime growth rate, Burgess et al. + (2026) Eq. (11). The paper gives the growth rate rather than an explicit width; + this is the corresponding length `delta_s_m` should sit within a few orders of magnitude of `d_beta` for a well-posed surface (`dels_db` is O(1)); a large gap flags a normalisation @@ -150,6 +170,32 @@ struct LayerWidths delta_s::ComplexF64 delta_s_m::Float64 d_beta::Float64 + delta_norm::Float64 + delta_visco::Float64 + delta_dr::Float64 +end + +""" + slayer_algebraic_widths(p::SLAYERParameters) -> (; delta_norm, delta_visco, delta_dr) + +The three comparison scales that follow from `p` alone, with no Riccati solve. `delta_dr` is +the channel the resistive-layer overlap criterion uses, so a scan needs only this. +""" +function slayer_algebraic_widths(p::SLAYERParameters) + # Layer normalization length: exactly 1/delta_n (delta_n = S^(1/3)/r_s), computed + # from rs and lu so it reads as a length. Burgess et al. (2026), the Δ̂' = Δ'/S^(1/3) + # normalization preceding Eq. (10). + delta_norm = p.rs * p.lu^(-1.0 / 3.0) + # Viscous-resistive broadening, P^(1/6) coefficient of Burgess et al. (2026) Eq. (11). + delta_visco = delta_norm * p.P_perp^(1.0 / 6.0) + # Diffusive-resistive width, Fitzpatrick (2025) Eq. (100), in meters. d_beta_hat is the + # normalized ion scale d_beta/rs; the shear is the r-based |s| SLAYER already carries. + # The paper's tau_A (its Eq. 74) carries no shear, whereas SLAYER's tau_h divides by n*s, + # so lu = tau_R/tau_h = (n|s|) * (tau_R/tau_A) supplies Eq. (100)'s explicit (n|s|)^(-1/2), + # leaving no separate shear factor in terms of lu. + d_beta_hat = p.d_beta / p.rs + delta_dr = d_beta_hat > 0 ? p.rs * p.lu^(-0.5) * p.P_perp^(0.25) / sqrt(d_beta_hat) : NaN + return (; delta_norm, delta_visco, delta_dr) end """ @@ -161,11 +207,16 @@ surface. Runs [`riccati_del_s`](@ref) for the dimensionless `δ_s / d_β` and scales by `p.d_beta` to obtain `δ_s` in meters. Keyword arguments are forwarded to `riccati_del_s`. + +The algebraic comparison scales `delta_norm`, `delta_visco` and `delta_dr` come from `p` +directly and do not depend on the Riccati solve; [`slayer_algebraic_widths`](@ref) returns +just those, for callers that do not need `delta_s`. """ function slayer_layer_thickness(p::SLAYERParameters; kwargs...) dels_db = riccati_del_s(p; kwargs...) delta_s = dels_db * p.d_beta + w = slayer_algebraic_widths(p) return LayerWidths(p.ising, p.m, p.n, dels_db, delta_s, abs(delta_s), - p.d_beta) + p.d_beta, w.delta_norm, w.delta_visco, w.delta_dr) end diff --git a/src/InnerLayer/SLAYER/SLAYER.jl b/src/InnerLayer/SLAYER/SLAYER.jl index 033f71eca..7f6c95c00 100644 --- a/src/InnerLayer/SLAYER/SLAYER.jl +++ b/src/InnerLayer/SLAYER/SLAYER.jl @@ -55,7 +55,7 @@ include("LayerInputs.jl") export SLAYERModel, SLAYERParameters, slayer_parameters export r_based_shear -export riccati_del_s, slayer_layer_thickness, LayerWidths +export riccati_del_s, slayer_layer_thickness, slayer_algebraic_widths, LayerWidths export surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs export NeoResistivityModel, SpitzerModel, SpitzerHarmModel, SauterNeoModel, RedlNeoModel diff --git a/src/Tearing/LayerInputs.jl b/src/Tearing/LayerInputs.jl index 3c1e2cc35..e7c5546ba 100644 --- a/src/Tearing/LayerInputs.jl +++ b/src/Tearing/LayerInputs.jl @@ -62,7 +62,7 @@ Throws if any surface's `restype` is still `nothing` — call `ForceFreeStates.resist_eval_all!(intr, equil)` first. """ function build_ggj_inputs(equil, sings, profiles::KineticProfiles; - mu_i::Real=2.0, zeff::Real=1.0, + mu_i::Real=2.0, zeff=1.0, v1_scale::Real=1.0, resistivity_model::NeoResistivityModel=SpitzerModel(), lnLambda_form::Symbol=:nrl) @@ -93,12 +93,13 @@ function build_ggj_inputs(equil, sings, profiles::KineticProfiles; # Shared Coulomb log and resistivity closure (identical to SLAYER # when the same resistivity_model is selected). lnLamb = coulomb_log_e(n_e, t_e; form=lnLambda_form) + zeff_k = zeff isa Real ? Float64(zeff) : Float64(zeff(sing.psifac)) if resistivity_model isa SpitzerModel - eta_use = eta_spitzer(n_e, t_e, zeff; lnLamb=lnLamb) + eta_use = eta_spitzer(n_e, t_e, zeff_k; lnLamb=lnLamb) else nuestar = nu_star_e(n_e, t_e, rg.R_major, rg.eps_local, - sing.q, zeff; lnLamb=lnLamb) - eta_use = eta_neoclassical(resistivity_model, n_e, t_e, zeff, + sing.q, zeff_k; lnLamb=lnLamb) + eta_use = eta_neoclassical(resistivity_model, n_e, t_e, zeff_k, rg.f_trap, nuestar; lnLamb=lnLamb) end rho = mu_i * M_P * n_e diff --git a/src/Tearing/LayerOverlap.jl b/src/Tearing/LayerOverlap.jl new file mode 100644 index 000000000..64c187d2d --- /dev/null +++ b/src/Tearing/LayerOverlap.jl @@ -0,0 +1,438 @@ +# LayerOverlap.jl +# +# Resistive-layer overlap scan: where adjacent rational surfaces' resistive layers run into each +# other. Past that point neither surface keeps a well-separated inner region, so the +# matched-asymptotic treatment is not defined there. Criterion: Fitzpatrick, Nucl. Fusion 2025 +# (doi 10.1088/1741-4326/ae4fdd) Sect. 5.9, judged with the diffusive-resistive width of its +# Eq. (100). How the result is used is documented in the integration-domain truncation section of +# docs/src/stability.md. +# +# Lives in `Tearing` rather than `InnerLayer/SLAYER` because it needs +# `ForceFreeStates._find_rational_surfaces`, and `InnerLayer` loads before `ForceFreeStates`. +# +# Surfaces beyond the equilibrium's own ψ grid are located on `Equilibrium.separatrix_q_law`: the +# logarithmic divergence of q at the X-point(s) on the last closed flux surface, with its +# location and slope taken from the 2D flux and its constant from the traced q on the outermost +# grid surface. The scan only recommends a domain; nothing here re-forms the equilibrium. +using ..Utilities: KineticProfiles +using ..ForceFreeStates: _find_rational_surfaces, SingType +using ..InnerLayer.SLAYER: SLAYERParameters, build_slayer_inputs, slayer_algebraic_widths, riccati_del_s, + surface_minor_radius, surface_da_dpsi, radial_label +using ..Utilities.NeoclassicalResistivity: NeoResistivityModel, SauterNeoModel +using ..Equilibrium: separatrix_q_law, separatrix_q, separatrix_q_psi, separatrix_q_dqdpsi, InverseIngest +using FastInterpolations: cubic_interp, ExtendExtrap, DerivOp +using Roots: find_zero, Brent +using Printf: @sprintf + +""" + LayerOverlapScan + +Per-surface resistive layer widths and the `psihigh` they imply, from +[`resistive_layer_overlap`](@ref). + +All widths are in **normalized flux**, converted from the metre-valued +`LayerWidths` scales by `Δψ = w / |da/dψ|` so they can be compared against +surface *separations* in ψ. Comparing the metre values directly against ψ is a +dimensional error. + +# Fields + + - `m`, `n` -- resonant mode numbers per surface, ordered by increasing ψ + - `psi` -- surface location in normalized flux + - `rs` -- minor radius at the surface in meters + - `delta_s_m` -- `|δ_s|`, the Riccati resistive layer width, in meters + - `width_delta_s` -- the same width in normalized flux, `|δ_s|/|da/dψ|` + - `width_visco` -- the visco-resistive comparison scale as a Δψ width + - `width_dr` -- the diffusive-resistive width of Fitzpatrick (2025) Eq. (100) as a Δψ + width. This is the criterion the recommendation uses: the paper's Sect. 5.6 shows that + strong shear near the separatrix forces every layer there into the DR regime + - `extrapolated` -- true when the surface was located outside the + equilibrium's ψ grid via the edge q-law + - `psihigh_delta_s`, `psihigh_visco`, `psihigh_dr` -- domain implied by each width + channel, `nothing` when that channel never overlaps within the scanned range + - `psihigh` -- the recommendation, equal to `psihigh_dr`; the other two channels are reported + for comparison and do not set the domain + - `first_overlap` -- index of the first surface that overlaps its inner + neighbour, `nothing` when none do + - `caveats` -- reasons the verdict itself is qualified: an overlap that could not be + established or excluded, or one decided by a surface placed on the separatrix law + - `notes` -- per-surface and bookkeeping detail that does not qualify the verdict +""" +struct LayerOverlapScan + m::Vector{Int} + n::Vector{Int} + psi::Vector{Float64} + rs::Vector{Float64} + delta_s_m::Vector{Float64} + width_delta_s::Vector{Float64} + width_visco::Vector{Float64} + width_dr::Vector{Float64} + extrapolated::Vector{Bool} + psihigh_delta_s::Union{Float64,Nothing} + psihigh_visco::Union{Float64,Nothing} + psihigh_dr::Union{Float64,Nothing} + psihigh::Union{Float64,Nothing} + first_overlap::Union{Int,Nothing} + caveats::Vector{String} + notes::Vector{String} +end + +# Keyword constructor, so an empty scan names only the fields it sets. +LayerOverlapScan(; m=Int[], n=Int[], psi=Float64[], rs=Float64[], delta_s_m=Float64[], + width_delta_s=Float64[], width_visco=Float64[], width_dr=Float64[], extrapolated=Bool[], + psihigh_delta_s=nothing, psihigh_visco=nothing, psihigh_dr=nothing, psihigh=nothing, + first_overlap=nothing, caveats=String[], notes=String[]) = + LayerOverlapScan(m, n, psi, rs, delta_s_m, width_delta_s, width_visco, width_dr, extrapolated, + psihigh_delta_s, psihigh_visco, psihigh_dr, psihigh, first_overlap, caveats, notes) + +# An inverse equilibrium is solved on a prescribed boundary and carries its own q out to ψ = 1, +# so the scan reads the real q there instead of extrapolating, and searches to ψ = 1 rather than +# stopping short of a separatrix it does not have. The direct/EFIT path cannot: GPEC recomputes q +# by field-line tracing only out to psihigh. +_inverse_q_spline(ingest::InverseIngest) = + cubic_interp(collect(Float64, ingest.sq_xs), collect(Float64, ingest.sq_fs[:, 3]); + extrap=ExtendExtrap()) + +# Locate q = q_target on the real q profile, between `psi_lo` and `psi_hi`. Returns nothing when +# the target is outside the range this equilibrium actually reaches. +function _real_q_surface(qspl, q_target::Real, psi_lo::Float64, psi_hi::Float64) + q_lo, q_hi = Float64(qspl(psi_lo)), Float64(qspl(psi_hi)) + (q_lo - q_target) * (q_hi - q_target) <= 0 || return nothing + return find_zero(p -> Float64(qspl(p)) - q_target, (psi_lo, psi_hi), Brent()) +end + +""" + resistive_layer_overlap(equil, profiles::KineticProfiles; n_tor, kwargs...) + -> LayerOverlapScan + +Locate the q = m/`n_tor` rational surfaces, compute each one's resistive layer +width, and report the outermost `psihigh` at which adjacent layers are still +separated. + +Surfaces inside the equilibrium grid come from +`ForceFreeStates._find_rational_surfaces` (Brent root-finding segmented between +q-extrema, so reverse shear is handled). Surfaces beyond the grid come from +[`separatrix_q_law`](@ref) when `extrapolate=true`, which is what lets the scan +recommend a domain *larger* than the one it was handed. That law diverges at the +X-point flux, which a reconstruction may place slightly outside its own boundary, +so rationals can lie at ψ > 1: inside the real separatrix but beyond the +equilibrium GPEC holds. Only the first of them is placed, to judge whether the +outermost surface inside is clean, and it is scored only if the kinetic profiles +extend that far. The recommendation never exceeds ψ = 1, and when no overlap can +be established the domain keeps its requested limit. + +Layer widths come from [`slayer_algebraic_widths`](@ref) via +[`build_slayer_inputs`](@ref), so the scan uses the same plasma inputs and the +same resistivity closure as the SLAYER analysis itself. The criterion channel needs no +Riccati solve; `delta_s` is solved per surface as a reported comparison only. + +# Arguments + + - `equil` -- a formed `PlasmaEquilibrium` + - `profiles` -- `KineticProfiles` spanning the scanned ψ range + +# Keyword arguments + + - `n_tor` -- toroidal mode number; surfaces are q = m/`n_tor` + + - `m_max` -- runaway backstop on the poloidal mode number (default `2000`), **not** the + physics limit: the outward search should terminate within `psi_gap` of the separatrix + flux, and if `m_max` is what stops it instead, the scan says so in `caveats` + because "no overlap found" would then rest on a non-physical cutoff + - `max_layer_solves` -- bound on how many surfaces get a Riccati layer solve (default `64`, + innermost kept). The first overlap decides the answer, so surfaces far beyond it only cost + time; when the bound bites, `caveats` records that a "no overlap" verdict is inconclusive + - `psihigh_safe` -- treat this as the outer edge of trusted equilibrium data + (default `equil.profiles.xs[end]`) + - `psi_gap` -- never place a surface closer than this to the separatrix flux on the + **direct** path (default `1e-9`, purely numerical: it stops the outward search from + crowding indefinitely onto the separatrix where surfaces become indistinguishable). + Inverse equilibria ignore it: they have a prescribed boundary and the search runs to + ψ = 1 on the real q. Geometry needs no clamp of its own -- `radial_label` derivatives are + analytic and valid under extrapolation + - `extrapolate` -- search past `psihigh_safe` on the separatrix q-law (default `true`) + - `theta` -- poloidal angle for the minor-radius chord (default `0.0`) + - `rs_method` -- radial label the widths and their ψ-conversion Jacobian are taken in + (default `:midplane`; see `radial_label`). The caller wiring this scan into the domain + truncation passes `:flux`, the label Fitzpatrick Eq. (100) is anchored to + - `bt`, `mu_i`, `zeff`, `chi_perp`, `chi_tor`, `resistivity_model`, `lnLambda_form` + -- passed through to `build_slayer_inputs`. The defaults are placeholders: a caller + scanning a real run must pass that run's own values, as the domain-truncation wiring does + +The overlap criterion is that surface `k` overlaps its inner neighbour when +`ψ_k − w_k/2 < ψ_{k-1} + w_{k-1}/2`. When that happens **both** surfaces are +contaminated, so the recommendation cuts at the inner edge of surface `k-1`, +leaving `k-2` as the outermost trustworthy surface. +""" +function resistive_layer_overlap(equil, profiles::KineticProfiles; + n_tor::Integer, + m_max::Integer=2000, + max_layer_solves::Integer=64, + psihigh_safe::Real=Float64(equil.profiles.xs[end]), + psi_gap::Real=1e-9, + extrapolate::Bool=true, + theta::Real=0.0, + rs_method::Symbol=:midplane, + bt=nothing, + mu_i::Real=2.0, + zeff=1.0, + chi_perp=1.0, + chi_tor=1.0, + resistivity_model::NeoResistivityModel=SauterNeoModel(), + lnLambda_form::Symbol=:nrl) + + n_tor > 0 || throw(ArgumentError("resistive_layer_overlap: n_tor must be positive, got $n_tor")) + psihigh_safe = Float64(psihigh_safe) + caveats, notes = String[], String[] + + # In-grid surfaces for this n, ordered outward. + found = [(m=s.m, psi=s.psifac, extrap=false) + for s in _find_rational_surfaces(equil, Int(n_tor), Int(n_tor)) + if s.psifac <= psihigh_safe && s.m <= m_max] + sort!(found; by=s -> s.psi) + + # Surfaces past the grid. An inverse equilibrium already knows its real q out to ψ = 1, so + # it is read rather than modelled; only the direct path needs the edge law. + edge_fit = nothing + inv_q = (equil.ingest isa InverseIngest) ? _inverse_q_spline(equil.ingest) : nothing + if extrapolate && inv_q !== nothing + psi_top = 1.0 # a prescribed boundary, not a separatrix + m_start = isempty(found) ? 1 : maximum(s.m for s in found) + 1 + n_before = length(found) + reached_end = false + for m in m_start:Int(m_max) + psi_m = _real_q_surface(inv_q, m / n_tor, psihigh_safe, psi_top) + if psi_m === nothing + reached_end = true # q never reaches m/n inside this equilibrium + break + end + psi_m > psihigh_safe || continue + push!(found, (m=m, psi=psi_m, extrap=true)) + end + !reached_end && m_start <= m_max && + push!(caveats, + "outward search stopped at the m_max = $m_max backstop rather than at the boundary " * + "q; a \"no overlap\" result here is not conclusive -- raise m_max") + length(found) > n_before && + push!(notes, "inverse equilibrium: surfaces beyond ψ = $(@sprintf("%.6f", psihigh_safe)) " * + "read from the real q profile out to ψ = 1 (no edge-law extrapolation)") + elseif extrapolate + edge_fit = separatrix_q_law(equil; psi_max=psihigh_safe) + if edge_fit === nothing + push!( + notes, + "no separatrix q-law: no X-point was found from the outermost flux surface, " * + "or the law it gives does not reproduce the traced shear " * + "there, so no surfaces are placed beyond ψ = $(@sprintf("%.6f", psihigh_safe))" + ) + else + # Every q = m/n above the anchor lies beyond the grid, so the search starts at the + # first such m whatever the in-grid search returned: no surface can fall between + # the two locators. + m_start = floor(Int, n_tor * abs(edge_fit.q_anchor)) + 1 + m_grid = isempty(found) ? 0 : maximum(s.m for s in found) + n_before = length(found) + # Terminates within psi_gap of the separatrix flux -- surfaces indistinguishable from + # it -- not on m_max. Reaching m_max means the scan was cut short for a non-physical + # reason and any "no overlap" verdict from it is unsafe, so record that. + reached_separatrix = false + for m in m_start:Int(m_max) + psi_m = separatrix_q_psi(edge_fit, sign(edge_fit.q_anchor) * m / n_tor; psi_gap=psi_gap) + if psi_m === nothing + reached_separatrix = true + break + end + push!(found, (m=m, psi=psi_m, extrap=true)) + # One rational past the equilibrium's boundary is enough: it decides whether the + # outermost surface inside is clean. Pairs further out lie wholly outside the + # equilibrium and cannot constrain a domain that ends at or before ψ = 1. + if psi_m > 1.0 + reached_separatrix = true + break + end + end + if !reached_separatrix && m_start <= m_max + push!(caveats, + "outward search stopped at the m_max = $m_max backstop rather than at the " * + "separatrix; a \"no overlap\" result here is not conclusive -- raise m_max") + end + 0 < m_grid < m_start - 1 && + push!(caveats, + "q = $(m_grid + 1)/$n_tor through $(m_start - 1)/$n_tor lie inside the grid but were " * + "not located there; an overlap between them may be missed") + length(found) > n_before && last(found).psi > 1.0 && + push!(notes, + "q = $(last(found).m)/$n_tor lies beyond the equilibrium's own boundary (ψ = 1), inside its " * + "X-point flux ψ = $(@sprintf("%.6f", edge_fit.psi_x[1])); it is placed only to judge whether the " * + "outermost surface inside is clean, and needs kinetic profiles past ψ = 1 to be scored") + end + end + + if isempty(found) + push!(notes, "no q = m/$n_tor surfaces found") + return LayerOverlapScan(; caveats=caveats, notes=notes) + end + + # Beyond the grid the equilibrium's q spline is a cubic extrapolation, so the radial label + # takes q there from whichever source placed the surfaces. + q_beyond = inv_q !== nothing ? (ψ -> Float64(inv_q(ψ))) : + edge_fit !== nothing ? (ψ -> separatrix_q(edge_fit, ψ)) : nothing + # Widths come back in metres of whichever radial label `rs_method` selects, so the + # conversion into normalised flux must use that same label's Jacobian. + _, _dr_dpsi = radial_label(equil; rs_method=rs_method, theta=theta, q_beyond=q_beyond) + + ms, ns, psis, rss = Int[], Int[], Float64[], Float64[] + dels_m, w_dels, w_visc, w_dr, extraps = Float64[], Float64[], Float64[], Float64[], Bool[] + unscored = Float64[] + + # Bound the per-surface Riccati work, and say so when the bound bites: a "no overlap" + # verdict from a truncated list is not conclusive. + if length(found) > max_layer_solves + push!(caveats, + "surface list truncated from $(length(found)) to the innermost $max_layer_solves for " * + "the layer solve (max_layer_solves); a \"no overlap\" verdict here is not conclusive") + found = found[1:max_layer_solves] + end + + for s in found + q_val = s.m / n_tor + # dq/dψ from the equilibrium spline in-grid, from the edge law outside it. + q1_val = if !s.extrap + Float64(equil.profiles.q_deriv(s.psi)) + elseif inv_q !== nothing + Float64(inv_q(s.psi; deriv=DerivOp(1))) # real profile, not a model + else + separatrix_q_dqdpsi(edge_fit, s.psi) + end + da_dpsi = _dr_dpsi(s.psi) + if !isfinite(da_dpsi) || da_dpsi == 0.0 + push!(notes, "m=$(s.m): da/dψ = $da_dpsi at ψ=$(@sprintf("%.6f", s.psi)); cannot convert width to flux units") + push!(unscored, s.psi) + continue + end + + sing = SingType(; m=[s.m], n=[Int(n_tor)], psifac=s.psi, rho=sqrt(s.psi), q=q_val, q1=q1_val) + # Built one surface at a time so a single degenerate surface (e.g. ω_*e == ω_*i, + # which makes iota_e singular) is recorded and skipped rather than aborting the scan. + params = try + build_slayer_inputs(equil, [sing], profiles; rs_method=rs_method, q_beyond=q_beyond, + bt=bt, mu_i=mu_i, zeff=zeff, chi_perp=chi_perp, chi_tor=chi_tor, + dr_val=0.0, dc_type=:none, theta=theta, + resistivity_model=resistivity_model, lnLambda_form=lnLambda_form) + catch err + # Any single-surface failure is recorded and skipped, never fatal. A DomainError is + # the common one: profiles that stop short of the extrapolated surfaces, or a + # negative-shear surface reaching a sqrt/fractional power in the layer parameters. + err isa Union{ArgumentError,DomainError} || rethrow() + msg = err isa ArgumentError ? err.msg : sprint(showerror, err) + push!(notes, "m=$(s.m) at ψ=$(@sprintf("%.6f", s.psi)): skipped -- $(replace(msg, '\n' => ' '))") + push!(unscored, s.psi) + continue + end + + # The criterion needs only the algebraic DR width, so the Riccati solve never gates it: + # dropping a surface here would leave _first_overlap_limit comparing its neighbours across + # a gap twice the physical one, which hides an overlap and reports too wide a domain. + alg = slayer_algebraic_widths(params[1]) + if !isfinite(alg.delta_dr) + push!( + notes, + "m=$(s.m) at ψ=$(@sprintf("%.6f", s.psi)): Eq. (100) width is not finite, so this " * + "surface could not be scored" + ) + push!(unscored, s.psi) + continue + end + # delta_s is reported for comparison only. A failure leaves it NaN, which makes its own + # channel report no overlap, and leaves the DR criterion untouched. + delta_s_m = try + w = abs(riccati_del_s(params[1]) * params[1].d_beta) + isfinite(w) || push!(notes, "m=$(s.m) at ψ=$(@sprintf("%.6f", s.psi)): del_s Riccati did not converge; " * + "delta_s channel only, criterion unaffected") + w + catch err + err isa Union{ArgumentError,DomainError} || rethrow() + push!(notes, "m=$(s.m) at ψ=$(@sprintf("%.6f", s.psi)): del_s Riccati failed; delta_s channel only, " * + "criterion unaffected") + NaN + end + lw = (; delta_s_m=delta_s_m, delta_dr=alg.delta_dr, delta_visco=alg.delta_visco) + + # Metres → normalized flux, so widths and surface separations are comparable. + push!(ms, s.m) + push!(ns, Int(n_tor)) + push!(psis, s.psi) + push!(rss, params[1].rs) + push!(dels_m, lw.delta_s_m) + push!(w_dels, lw.delta_s_m / abs(da_dpsi)) + push!(w_visc, lw.delta_visco / abs(da_dpsi)) + push!(w_dr, lw.delta_dr / abs(da_dpsi)) + push!(extraps, s.extrap) + end + + ph_dels, _ = _first_overlap_limit(psis, w_dels) + ph_visc, _ = _first_overlap_limit(psis, w_visc) + ph_dr, k_dr = _first_overlap_limit(psis, w_dr) + # The DR channel is the criterion (see the `width_dr` field above); the other two are + # reported for comparison only. + recommended = ph_dr + # A surface that could not be scored hides any overlap with its neighbours, so it matters + # wherever it sits inside the reported overlap point, or anywhere when none was found. + n_hidden = count(p -> k_dr === nothing || p < psis[k_dr], unscored) + n_hidden > 0 && + push!(caveats, + k_dr === nothing ? + "no overlap was established, with $n_hidden surface(s) that could not be scored (typically " * + "kinetic profiles that stop short of them); the domain keeps its requested limit" : + "$n_hidden surface(s) inside the reported overlap point could not be scored (typically " * + "kinetic profiles that stop short of them); the first overlap may lie further in, so " * + "this result is an outer bound") + k_dr !== nothing && (extraps[k_dr] || extraps[k_dr-1]) && + push!(caveats, + "the overlap point is decided by a surface placed on the separatrix law (q = $(ms[k_dr-1])/$n_tor " * + "against $(ms[k_dr])/$n_tor), not by one traced inside the grid") + + return LayerOverlapScan(ms, ns, psis, rss, dels_m, w_dels, w_visc, w_dr, extraps, + ph_dels, ph_visc, ph_dr, recommended, k_dr, caveats, notes) +end + +# Walk outward and return (recommended psilim, index of first overlapping surface). +# Overlap at k contaminates both k and k-1, so the limit is surface k-1's own position +# (Fitzpatrick 2025 Sect. 5.9 retains 0 < psi < 1 - eps_c with psi[k-1] as 1 - eps_c); +# offsetting by ±w[k-1]/2 either drops the last clean surface or can exceed psi = 1. +function _first_overlap_limit(psi::Vector{Float64}, w::Vector{Float64}) + for k in 2:length(psi) + if psi[k] - w[k] / 2 < psi[k-1] + w[k-1] / 2 + return (psi[k-1], k) + end + end + return (nothing, nothing) +end + +function Base.show(io::IO, ::MIME"text/plain", s::LayerOverlapScan) + println(io, "LayerOverlapScan: $(length(s.psi)) surface(s), n = ", + isempty(s.n) ? "-" : string(s.n[1])) + if !isempty(s.psi) + println(io, rpad("q", 8), rpad("psi", 13), rpad("r_s [m]", 11), + rpad("dpsi(del_s)", 14), rpad("dpsi(visco)", 14), rpad("dpsi(DR)", 14), "source") + for i in eachindex(s.psi) + println(io, + rpad("$(s.m[i])/$(s.n[i])", 8), + rpad(@sprintf("%.7f", s.psi[i]), 13), + rpad(@sprintf("%.5f", s.rs[i]), 11), + rpad(@sprintf("%.4e", s.width_delta_s[i]), 14), + rpad(@sprintf("%.4e", s.width_visco[i]), 14), + rpad(@sprintf("%.4e", s.width_dr[i]), 14), + s.extrapolated[i] ? "extrapolated" : "in-grid") + end + end + _fmt(x) = x === nothing ? "none (no overlap in range)" : @sprintf("%.7f", x) + println(io, " psihigh from |del_s| : ", _fmt(s.psihigh_delta_s)) + println(io, " psihigh from visco : ", _fmt(s.psihigh_visco)) + println(io, " psihigh from Eq.(100): ", _fmt(s.psihigh_dr)) + println(io, " recommended psihigh : ", _fmt(s.psihigh)) + for note in s.notes + println(io, " note: ", note) + end + return nothing +end diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 1a12ff7e6..41b862a5d 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -23,7 +23,7 @@ constructor. (multi-surface determinant) - `dc_type` -- critical-Δ offset selector, one of `:none`, `:lar`, `:rfitzp`, `:toroidal` (χ_‖-matching critical-Δ formulas, - Connor-Hastie-Helander 2015) + Connor, Ham, Hastie & Liu 2015) - `msing_max` -- number of surfaces to include in the coupled determinant (default 3; capped at `length(sings)` at runtime) @@ -33,7 +33,13 @@ constructor. `B_T = F(ψ)/(2π·R₀)` per surface from the equilibrium's F-spline; a scalar or a callable of `psi` overrides it - `mu_i` -- ion mass in proton-mass units (default 2.0 for D) - - `zeff` -- effective charge + - `zeff` -- fallback effective charge. `Z_eff(ψ)` is taken, in order of precedence, from + the kinetic file's own `zeff` profile; else derived from the file's `n_i/n_e` with the + impurity charge `zimp`; else this scalar at every surface. A dataset that is absent or + all-zero counts as not provided + - `zimp` -- charge of the single impurity species that closes quasi-neutrality when + `Z_eff(ψ)` is derived from `n_i/n_e` (default 6.0, carbon); the main ion has charge 1. + Unused when the file supplies a `zeff` profile - `chi_perp`, `chi_tor` -- fallback perpendicular / toroidal heat diffusivity [m²/s], used only when the kinetic file carries no usable `chi_e`/`chi_phi` profile (dataset absent or all-zero); otherwise the @@ -107,6 +113,7 @@ there is one consistent interface for resistive and kinetic profiles. bt::Union{Float64,Nothing} = nothing mu_i::Float64 = 2.0 zeff::Float64 = 1.0 + zimp::Float64 = 6.0 chi_perp::Float64 = 1.0 chi_tor::Float64 = 1.0 dr_val::Union{Float64,Nothing} = nothing @@ -183,6 +190,8 @@ function validate(ctrl::SLAYERControl) ctrl.lnLambda_form in _VALID_LNLAMBDA_FORMS || throw(ArgumentError("SLAYERControl: lnLambda_form=$(ctrl.lnLambda_form) " * "not in $(_VALID_LNLAMBDA_FORMS)")) + ctrl.zimp > 1 || + throw(ArgumentError("SLAYERControl: zimp=$(ctrl.zimp) must exceed the main-ion charge 1")) ctrl.msing_max >= 1 || throw(ArgumentError("SLAYERControl: msing_max=$(ctrl.msing_max) must be ≥ 1")) ctrl.nre >= 2 && ctrl.nim >= 2 || diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index b0fc8d2e4..87330a52e 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -90,10 +90,10 @@ const TEARING_H5_ANNOTATIONS = [ "PerSurface/sval_r" => (; long_name="r-based magnetic shear r_s·(dq/dr)/q (Fitzpatrick convention)", dims=("surface",)), "PerSurface/D_R" => (; long_name="resistive interchange D_R = E + F + H² for the critical-Δ formula (auto-derived from GGJ coefficients unless overridden)", dims=("surface",)), - "PerSurface/D_geo" => (; long_name="Connor-Hastie-Helander 2015 Eq. 59 geometric factor (0 unless supplied)", dims=("surface",)), + "PerSurface/D_geo" => (; long_name="Connor et al. 2015 Eq. 59 geometric factor (0 unless supplied)", dims=("surface",)), "PerSurface/eta" => (; long_name="parallel resistivity at each surface", units="Ohm*m", dims=("surface",)), "PerSurface/d_beta" => (; long_name="β-weighted ion drift scale d_β", units="m", dims=("surface",)), - "PerSurface/D_c_offset" => (; long_name="critical-Δ offset from χ_∥/χ_⊥ matching (Connor-Hastie-Helander 2015 Eq. 59)", dims=("surface",)), + "PerSurface/D_c_offset" => (; long_name="critical-Δ offset from χ_∥/χ_⊥ matching (Connor et al. 2015 Eq. 59)", dims=("surface",)), "PerSurface/D_c_type" => (; long_name="per-surface D_c prescription label", dims=("surface",)), "PerSurface/k_ref" => (; long_name="reference-length ratio K = r_s·(dψ_N/dr) at each surface", dims=("surface",)), "PerSurface/alpha_mercier" => (; long_name="Mercier Frobenius exponent α = √(−D_I) at each surface (Glasser-Greene-Johnson 1975 Eq. 48)", dims=("surface",)), diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 2459087cf..0b2758293 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -20,7 +20,10 @@ # reader and adapt them to the SLAYER inner layer. Returns the spline-based # `KineticProfiles` plus χ⊥(ψ)/χ_φ(ψ) callables built from the file's # `chi_e`/`chi_phi` (or `nothing` when the file carries no usable χ, in which -# case the caller falls back to the scalar `control.chi_perp`/`chi_tor`). +# case the caller falls back to the scalar `control.chi_perp`/`chi_tor`), and +# a Z_eff(ψ) callable: the file's own `zeff` profile when it carries one, else +# derived from its `n_i/n_e` (or `nothing` when the file carries neither, in which +# case the caller falls back to the scalar `control.zeff`). function _load_profiles(control::SLAYERControl, dir_path::AbstractString) isempty(control.profile_file) && error("run_slayer: [SLAYER] profile_file is empty — point it at a " * @@ -53,7 +56,25 @@ function _load_profiles(control::SLAYERControl, dir_path::AbstractString) cubic_interp(psi_xs, collect(Float64, v)) chi_perp = _chi_spline(data.chi_e) chi_tor = _chi_spline(data.chi_phi) - return (profiles=profiles, chi_perp=chi_perp, chi_tor=chi_tor) + # Z_eff(ψ): a profile supplied in the file wins, since a real impurity mix is not recoverable + # from n_i/n_e. Otherwise quasi-neutrality with one impurity species of charge `zimp` and a + # charge-1 main ion, the single-ion closure the kinetic spline set uses: + # Z_eff = zimp − (n_i/n_e)(zimp − 1). + _usable(v) = v !== nothing && !all(iszero, v) + zeff = if _usable(data.zeff) + minimum(data.zeff) >= 1.0 - 1e-8 || + error("run_slayer: kinetic file '$path' has zeff = $(minimum(data.zeff)) < 1 — unphysical") + cubic_interp(psi_xs, collect(Float64, data.zeff)) + elseif !_usable(data.n_i) + nothing + else + z = [ne > 0 ? control.zimp - (ni / ne) * (control.zimp - 1.0) : control.zimp for (ni, ne) in zip(data.n_i, data.n_e)] + minimum(z) >= 1.0 - 1e-8 || + error("run_slayer: kinetic file '$path' has n_i > n_e, giving Z_eff = $(minimum(z)) < 1 — unphysical") + cubic_interp(psi_xs, z) + end + zeff_source = _usable(data.zeff) ? :file : zeff === nothing ? :scalar : :closure + return (profiles=profiles, chi_perp=chi_perp, chi_tor=chi_tor, zeff=zeff, zeff_source=zeff_source) end # --------------------------------------------------------------------- @@ -437,13 +458,22 @@ function run_slayer(equil, surfaces::AbstractVector, delta_prime_matrix::Abstrac loaded = _load_profiles(control, dir_path) profiles = loaded.profiles + # Z_eff(ψ) from the kinetic file when it carries n_i, else the scalar fallback. + zeff = loaded.zeff === nothing ? control.zeff : loaded.zeff + if loaded.zeff_source === :scalar + @warn "SLAYER: kinetic file has no usable zeff or n_i profile (dataset absent or all-zero); " * + "using the scalar control.zeff = $(control.zeff) at every surface." + else + zeff_origin = loaded.zeff_source === :file ? "zeff profile" : "n_i/n_e with zimp = $(control.zimp)" + @info "SLAYER: Z_eff(ψ) from the kinetic file's $zeff_origin" + end if control.inner_model in (:ggj_shooting, :ggj_galerkin) # GGJ γ-extraction is future work; `run_slayer_from_inputs` emits the # warning once the model is built (so direct callers see it too). params = build_ggj_inputs(equil, surfaces, profiles; mu_i=control.mu_i, - zeff=control.zeff, + zeff=zeff, resistivity_model=_build_resistivity_model(control.resistivity_model), lnLambda_form=control.lnLambda_form) else @@ -463,7 +493,7 @@ function run_slayer(equil, surfaces::AbstractVector, delta_prime_matrix::Abstrac params = build_slayer_inputs(equil, surfaces, profiles; bt=bt, mu_i=control.mu_i, - zeff=control.zeff, + zeff=zeff, chi_perp=chi_perp, chi_tor=chi_tor, dr_val=control.dr_val, diff --git a/src/Tearing/Tearing.jl b/src/Tearing/Tearing.jl index 745a30857..80bce7b23 100644 --- a/src/Tearing/Tearing.jl +++ b/src/Tearing/Tearing.jl @@ -22,6 +22,7 @@ using ..Utilities import ..InnerLayer as InnerLayer include("LayerInputs.jl") +include("LayerOverlap.jl") include("Dispersion/Dispersion.jl") include("Runner/Runner.jl") @@ -30,5 +31,6 @@ import .Runner as Runner export InnerLayer, Dispersion, Runner export build_ggj_inputs +export resistive_layer_overlap, LayerOverlapScan end # module Tearing diff --git a/test/runtests.jl b/test/runtests.jl index e1256e56d..e18762e18 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -49,6 +49,7 @@ else include("./runtests_slayer_params.jl") include("./runtests_slayer_riccati.jl") include("./runtests_slayer_inputs.jl") + include("./runtests_layer_overlap_fitzpatrick.jl") include("./runtests_dispersion_residual.jl") include("./runtests_dispersion_coupled.jl") include("./runtests_dispersion_coupled_full.jl") diff --git a/test/runtests_grid_refinement.jl b/test/runtests_grid_refinement.jl index 35c350c90..66d34332a 100644 --- a/test/runtests_grid_refinement.jl +++ b/test/runtests_grid_refinement.jl @@ -182,5 +182,26 @@ const GridRef = GeneralizedPerturbedEquilibrium.Equilibrium equil2 = GridRef.setup_equilibrium(eq_config, sol_config; override_psi_nodes=grid) @test equil2.profiles.xs == grid @test equil2.rzphi_xs == grid + + # Solovev has finite edge q, so the auto grid must classify it as limited: no separatrix + # log law, on either pass, even with the absolute r² floor disabled (the relative + # log-vs-linear comparison is what rejects it). + @test GridRef.edge_q_law(equil1) === nothing + @test GridRef.edge_q_law(equil2) === nothing + @test GridRef.edge_q_law(equil1; min_r2=-Inf) === nothing + end + + @testset "edge_q_law accepts a diverted edge" begin + using TOML + dir_path = joinpath(dirname(@__DIR__), "examples", "DIIID-like_ideal_example") + inputs = TOML.parsefile(joinpath(dir_path, "gpec.toml")) + equil = GridRef.setup_equilibrium(GridRef.EquilibriumConfig(inputs["Equilibrium"], dir_path), nothing) + fit = GridRef.edge_q_law(equil) + @test fit !== nothing + @test fit.A < 0 + @test fit.r2_log > fit.r2_linear + # The shear-consistency check is live: this fit's slope is ~5% off the measured + # -(1-psi)q', so a tolerance tighter than that must reject it. + @test GridRef.edge_q_law(equil; shear_ratio=1.01) === nothing end end diff --git a/test/runtests_layer_overlap_fitzpatrick.jl b/test/runtests_layer_overlap_fitzpatrick.jl new file mode 100644 index 000000000..64bfffed3 --- /dev/null +++ b/test/runtests_layer_overlap_fitzpatrick.jl @@ -0,0 +1,351 @@ +# External-reference check of the diffusive-resistive layer width against a published +# result: Fitzpatrick, Nucl. Fusion 2025 (doi 10.1088/1741-4326/ae4fdd), Sect. 5.8-5.9. +# +# The paper reports that for its model JET equilibrium the resistive layers of adjacent +# rational surfaces first overlap at Psi = 0.9985 for n = 1 (its Fig. 9) and Psi = 0.9952 +# for n = 4 (its Fig. 10). This rebuilds that equilibrium from the paper's own definitions +# and drives GPEC's `delta_dr` through `slayer_parameters`, so the width formula and the +# tau_R / tau_A / d_beta chain feeding it are checked against a number in print rather than +# against GPEC's own history. +@testset "Layer overlap vs Fitzpatrick (2025) JET model" begin + using GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER: slayer_parameters, + slayer_layer_thickness, slayer_algebraic_widths, SpitzerHarmModel + using QuadGK: quadgk + using Roots: find_zero, Bisection + + # Sect. 5.8 parameters. + B0, R0, A_MIN = 3.45, 2.96, 1.25 + ZEFF, MNUM, CHI = 10.0, 2.0, 1.0 + Q0, Q95, Q105 = 1.01, 3.5, 4.0 + + # Model safety factor, Eqs. (32)-(35). alpha- depends on rhat95, which depends on the + # profile through Eq. (36), so the pair is solved self-consistently. + q_in(r, am) = Q0 - am * log(1 - r^2) + q_out(r, ap) = -ap * log(r^2 - 1) + q_any(r, am, ap) = r < 1 ? q_in(r, am) : q_out(r, ap) + psi_raw(r, am, ap) = quadgk(x -> x / q_any(x, am, ap), 0.0, r; rtol=1e-11)[1] + + am, ap = 1.0, 1.0 + fp_converged = false + for _ in 1:200 + psi_sep = psi_raw(1.0 - 1e-12, am, ap) + r95 = find_zero(r -> psi_raw(r, am, ap) / psi_sep - 0.95, (0.5, 1 - 1e-9), Bisection()) + f105(r) = (psi_sep + quadgk(x -> x / q_out(x, ap), 1 + 1e-12, r; rtol=1e-10)[1]) / psi_sep - 1.05 + hi = 1.0 + 1e-9 + while hi < 1.40 && f105(hi) < 0 # q_out turns negative past r = sqrt(2) + hi += 0.005 + end + hi < 1.40 || error("Fitzpatrick JET model: failed to bracket the Psi = 1.05 surface before q_out degenerates") + r105 = find_zero(f105, (1 + 1e-9, hi), Bisection()) + am_new = -(Q95 - Q0) / log(1 - r95^2) + ap_new = -Q105 / log(r105^2 - 1) + fp_converged = abs(am_new - am) < 1e-12 && abs(ap_new - ap) < 1e-12 + am, ap = am_new, ap_new + fp_converged && break + end + fp_converged || error("Fitzpatrick JET model: the alpha-/alpha+ fixed point did not converge in 200 iterations") + psi_sep = psi_raw(1.0 - 1e-12, am, ap) + q_of(r) = q_any(r, am, ap) + Psi_of(r) = r < 1 ? psi_raw(r, am, ap) / psi_sep : + (psi_sep + quadgk(x -> x / q_out(x, ap), 1 + 1e-12, r; rtol=1e-10)[1]) / psi_sep + dq_dr(r) = r < 1 ? am * 2r / (1 - r^2) : -ap * 2r / (r^2 - 1) + s_of(r) = r * dq_dr(r) / q_of(r) + + # mtanh edge profiles read off the paper's Fig. 8, anchored at Psi = 0.96 and at the + # separatrix (Sect. 5.9 states Te ~ 100 eV there). + function make_tanh(f096, f100, f_ped, f_sol) + a = (f_ped - f_sol) / 2 + x1 = atanh(clamp(1 - (f096 - f_sol) / a, -0.999, 0.999)) + x2 = atanh(clamp(1 - (f100 - f_sol) / a, -0.999, 0.999)) + d = 0.04 / (x2 - x1) + p0 = 0.96 - x1 * d + return P -> f_sol + a * (1 - tanh((P - p0) / d)) + end + te_of = make_tanh(400.0, 100.0, 600.0, 20.0) + ne_of = make_tanh(3.5e19, 1.2e19, 5.0e19, 3.0e18) + + # Width in units of rhat, through the shipped code path. Eq. (100) is algebraic, so this takes + # the no-ODE path; the two are asserted equal once below. + function params_at(r, n) + P = Psi_of(r) + te = te_of(P) + ne = ne_of(P) + p = slayer_parameters(; n_e=ne, t_e=te, t_i=te, + omega=0.0, omega_e=1.0e4, omega_i=5.0e3, + qval=q_of(r), sval_r=s_of(r), bt=B0, + rs=r * A_MIN, R0=R0, mu_i=MNUM, zeff=ZEFF, + chi_perp=CHI, chi_tor=CHI, m=1, n=n, + resistivity_model=SpitzerHarmModel(), lnLambda_form=:nrl) + return p + end + width_hat(r, n) = slayer_algebraic_widths(params_at(r, n)).delta_dr / A_MIN + + # The scan's criterion channel must be exactly what the full width routine reports, so the + # no-ODE path cannot drift from it. + let p = params_at(0.9, 1) + @test slayer_algebraic_widths(p).delta_dr == slayer_layer_thickness(p).delta_dr + @test slayer_algebraic_widths(p).delta_visco == slayer_layer_thickness(p).delta_visco + end + + function surfaces(n; rmin=0.90, rmax=1.35) + out = Tuple{Int,Float64}[] + for m in 1:400 + qt = m / n + qt < q_of(rmin) && continue + r = if qt < q_of(1 - 1e-13) + find_zero(x -> q_of(x) - qt, (rmin, 1 - 1e-13), Bisection()) + elseif qt > q_of(1 + 1e-13) && qt < q_of(rmax) + find_zero(x -> q_of(x) - qt, (1 + 1e-13, rmax), Bisection()) + else + nothing + end + r === nothing || push!(out, (m, r)) + end + sort!(out; by=t -> t[2]) + return out + end + + # Inner boundary of the overlap region: the first surface whose layer runs into its + # inner neighbour, reported at that neighbour's location. + function overlap_psi(n) + S = surfaces(n) + for k in 2:length(S) + gap = S[k][2] - S[k-1][2] + (width_hat(S[k-1][2], n) + width_hat(S[k][2], n)) / 2 >= gap && return Psi_of(S[k-1][2]) + end + return nothing + end + + # The self-consistent profile must reproduce the paper's own inputs before the widths + # mean anything. + @test isapprox(q_of(0.0), Q0; atol=1e-10) + @test am > 0 && ap > 0 + + psi1 = overlap_psi(1) + @test psi1 !== nothing + psi1 === nothing || @info "Fitzpatrick 2025 overlap, n = 1" psi=psi1 paper=0.9985 deviation=abs(psi1 - 0.9985) atol=1e-4 + # Paper: Psi = 0.9985 (n = 1); measured deviation 1.1e-5. The bound covers the pedestal read + # off Fig. 8 and the resistivity closure (Spitzer-Harm here against the paper's Eqs. 70-72), + # and still leaves ~9x margin. Pure quadrature and root-finding, no BLAS, so platform spread + # sits far below this. + @test isapprox(psi1, 0.9985; atol=1e-4) + + psi4 = overlap_psi(4) + @test psi4 !== nothing + psi4 === nothing || @info "Fitzpatrick 2025 overlap, n = 4" psi=psi4 paper=0.9952 deviation=abs(psi4 - 0.9952) atol=1e-3 + # Paper: Psi = 0.9952 (n = 4); measured deviation 2.3e-4. Looser than n = 1 because these + # surfaces sit further in, where the digitised pedestal shape matters more. + @test isapprox(psi4, 0.9952; atol=1e-3) + # The overlap region must move inward with n, which is the paper's reported scaling. + @test psi4 < psi1 +end + +# Wiring of the overlap cap into the integration domain. The scan itself is checked against the +# paper above; this checks that the cap reaches `sing_lim!` and behaves as an upper bound. +@testset "Layer-overlap cap wiring" begin + using GeneralizedPerturbedEquilibrium.ForceFreeStates: sing_lim!, ForceFreeStatesInternal, + ForceFreeStatesControl + using GeneralizedPerturbedEquilibrium + using GeneralizedPerturbedEquilibrium.Equilibrium + using TOML + + dir_path = joinpath(dirname(@__DIR__), "examples", "Solovev_ideal_example") + inputs = TOML.parsefile(joinpath(dir_path, "gpec.toml")) + equil = Equilibrium.setup_equilibrium( + Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dir_path), + Equilibrium.SolovevConfig(inputs["SOL_INPUT"])) + + _fresh() = (i = ForceFreeStatesInternal(); i.nlow = 1; i.nhigh = 1; i) + ctrl = ForceFreeStatesControl(; set_psilim_via_dmlim=false, qhigh=1e3) + + # Baseline: no cap. + base = _fresh() + sing_lim!(base, ctrl, equil) + + # A cap inside the domain must pull qlim down and move psilim inward. + inside = _fresh() + cap = 0.5 * (equil.profiles.xs[1] + base.psilim) + sing_lim!(inside, ctrl, equil; psilim_cap=cap) + @test inside.qlim < base.qlim + @test inside.psilim < base.psilim + @test isapprox(inside.psilim, cap; rtol=1e-6) + + # A cap beyond psihigh is inert by construction -- the bound never widens the domain. + beyond = _fresh() + sing_lim!(beyond, ctrl, equil; psilim_cap=base.psilim + 1e-3) + @test beyond.qlim == base.qlim + @test beyond.psilim == base.psilim + + # nothing is the same as not passing it at all. + none = _fresh() + sing_lim!(none, ctrl, equil; psilim_cap=nothing) + @test none.qlim == base.qlim + @test none.psilim == base.psilim + + # With dmlim on, the (last rational + dmlim)/n offset must never step past qhigh or the cap. + ctrl_dm = ForceFreeStatesControl(; set_psilim_via_dmlim=true, dmlim=0.2, qhigh=1e3) + dm_base = _fresh() + sing_lim!(dm_base, ctrl_dm, equil) + N = floor(Int, equil.params.qmax - 0.15) + @test N >= 2 + for qh in (Float64(N), N + 0.1) + c = ForceFreeStatesControl(; set_psilim_via_dmlim=true, dmlim=0.2, qhigh=qh) + i = _fresh() + sing_lim!(i, c, equil) + @test i.qlim ≈ N - 1 + 0.2 + @test i.qlim <= qh + end + + # A cap placed exactly on a rational must exclude it on either side of the float knife-edge. + q_of = equil.profiles.q_spline + lo, hi = equil.profiles.xs[1], equil.params.psihigh_resolved + for _ in 1:200 + mid = 0.5 * (lo + hi) + q_of(mid) < N ? (lo = mid) : (hi = mid) + end + psi_N = 0.5 * (lo + hi) + for cap_k in (prevfloat(psi_N, 1000), psi_N, nextfloat(psi_N, 1000)) + i = _fresh() + sing_lim!(i, ctrl_dm, equil; psilim_cap=cap_k) + @test i.qlim ≈ N - 1 + 0.2 + @test i.psilim <= cap_k + end + @test dm_base.qlim > N - 1 + 0.2 + + # The control flag exists and is opt-in. + @test ForceFreeStatesControl().psilim_from_layer_overlap == false + + # Multi-n: the innermost overlap point binds; scans with no overlap never do. + Scan = GeneralizedPerturbedEquilibrium.Tearing.LayerOverlapScan + mkscan(ph) = Scan(; psihigh_dr=ph, psihigh=ph, first_overlap=(ph === nothing ? nothing : 1)) + binding = GeneralizedPerturbedEquilibrium._binding_overlap + s1, s2, s_none = mkscan(0.9997), mkscan(0.9952), mkscan(nothing) + @test binding([s1, s2, s_none]) === s2 + @test binding([s_none, s1]) === s1 + @test binding([s_none, mkscan(nothing)]) === s_none + @test binding([s1]) === s1 + + # The scan must see the run's own plasma: every [SLAYER] setting the layer width depends on + # reaches it, and changing one changes the widths. + trunc_dir = joinpath(dirname(@__DIR__), "examples", "DIIID-like_truncation_example") + slayer_section = TOML.parsefile(joinpath(trunc_dir, "gpec.toml"))["SLAYER"] + slayer_inputs(over) = Dict{String,Any}("SLAYER" => merge(slayer_section, over)) + overlap_inputs(over; dir=trunc_dir) = GeneralizedPerturbedEquilibrium._overlap_slayer_inputs(slayer_inputs(over), dir) + oi = overlap_inputs(Dict("mu_i" => 3.0, "bt" => 1.7, "theta_sample" => 0.3, "lnLambda_form" => "sauter")) + @test (oi.mu_i, oi.bt, oi.theta, oi.lnLambda_form) == (3.0, 1.7, 0.3, :sauter) + + # Z_eff(psi) comes from the kinetic file's own n_i/n_e with the [SLAYER] impurity charge. + kin_path = joinpath(trunc_dir, slayer_section["profile_file"]) + kin = Equilibrium.read_kinetic_file(kin_path) + j = length(kin.psi) ÷ 2 + @test oi.zeff(kin.psi[j]) ≈ 6.0 - kin.n_i[j] / kin.n_e[j] * 5.0 + @test overlap_inputs(Dict("zimp" => 10.0)).zeff(kin.psi[j]) ≈ 10.0 - kin.n_i[j] / kin.n_e[j] * 9.0 + @test 1.0 < oi.zeff(kin.psi[j]) < 6.0 + # A zeff profile supplied in the file overrides the n_i/n_e closure; a file with neither falls + # back to the scalar [SLAYER] zeff. + mktempdir() do tmp + own = fill(3.2, length(kin.psi)) + base = (; psi=kin.psi, n_e=kin.n_e, T_i=kin.T_i, T_e=kin.T_e, omega_E=kin.omega_E, chi_e=kin.chi_e, chi_phi=kin.chi_phi) + Equilibrium.write_kinetic_h5(joinpath(tmp, "own_zeff.h5"), Equilibrium.KineticProfileData(; base..., n_i=kin.n_i, zeff=own)) + @test Equilibrium.read_kinetic_file(joinpath(tmp, "own_zeff.h5")).zeff == own + @test overlap_inputs(Dict("profile_file" => "own_zeff.h5", "zimp" => 10.0, "zeff" => 2.5); dir=tmp).zeff(kin.psi[j]) ≈ 3.2 + Equilibrium.write_kinetic_h5(joinpath(tmp, "no_ni.h5"), Equilibrium.KineticProfileData(; base...)) + @test overlap_inputs(Dict("profile_file" => "no_ni.h5", "zeff" => 2.5); dir=tmp).zeff == 2.5 + end + + widths(over) = GeneralizedPerturbedEquilibrium._layer_overlap_scan(overlap_inputs(over), _fresh(), equil).width_dr + w_ref = widths(Dict{String,Any}()) + @test !isempty(w_ref) && all(isfinite, w_ref) + for over in (Dict("zimp" => 10.0), Dict("resistivity_model" => "spitzer"), Dict("lnLambda_form" => "sauter")) + @test all(.!isapprox.(widths(over), w_ref; rtol=1e-3)) + end + # Ion mass and field strength cancel between the Lundquist number and the ion sound radius in + # Eq. (100), so the criterion width does not depend on them; the scalar zeff is only a + # fallback, so it is inert while the file carries n_i. + for over in (Dict("mu_i" => 1.0), Dict("bt" => 1.0), Dict("zeff" => 3.0)) + @test all(isapprox.(widths(over), w_ref; rtol=1e-10)) + end + + # The no-surfaces path must return an empty scan with a note, not throw (m_max=0 with + # extrapolation off forces it deterministically). + @testset "resistive_layer_overlap: empty surface list" begin + using GeneralizedPerturbedEquilibrium: Tearing, Utilities + npsi = 8 + psi_grid = collect(range(0.0, 1.0; length=npsi)) + profs = Utilities.KineticProfiles(; psi=psi_grid, + n_e=fill(1.0e19, npsi), T_e=fill(500.0, npsi), T_i=fill(500.0, npsi), + omega=zeros(npsi), omega_e=zeros(npsi), omega_i=zeros(npsi)) + scan = Tearing.resistive_layer_overlap(equil, profs; n_tor=1, m_max=0, extrapolate=false) + @test scan isa Tearing.LayerOverlapScan + @test isempty(scan.psi) + @test scan.psihigh === nothing && scan.first_overlap === nothing + @test any(contains("no q ="), scan.notes) + end +end + +# Surfaces beyond the grid are placed on the X-point law. Checked against the equilibrium's own +# traced q: a law built on a short grid must land where a longer grid actually finds the surfaces. +@testset "Separatrix q-law" begin + using GeneralizedPerturbedEquilibrium.Equilibrium: separatrix_q_law, separatrix_q, separatrix_q_psi, separatrix_q_dqdpsi + using GeneralizedPerturbedEquilibrium.ForceFreeStates: ForceFreeStatesInternal + + # The functions themselves, on a two-X-point law with known parameters. + law2 = (psi_x=[1.0002, 1.01], A=[0.5, 0.3], psi_anchor=0.995, q_anchor=4.3) + @test separatrix_q(law2, law2.psi_anchor) == law2.q_anchor + for q_target in (4.5, 6.0, 9.0) + psi_q = separatrix_q_psi(law2, q_target) + @test law2.psi_anchor < psi_q < law2.psi_x[1] + @test separatrix_q(law2, psi_q) ≈ q_target rtol = 1e-10 + end + @test separatrix_q_dqdpsi(law2, 0.999) ≈ 0.5 / (1.0002 - 0.999) + 0.3 / (1.01 - 0.999) + @test separatrix_q_psi(law2, 4.3) === nothing # not beyond the anchor + @test separatrix_q_psi(law2, 4.0) === nothing + @test separatrix_q_psi(law2, 40.0) === nothing # closer to the separatrix than psi_gap + @test separatrix_q_psi(law2, 40.0; psi_gap=1e-40) !== nothing + + trunc_dir = joinpath(dirname(@__DIR__), "examples", "DIIID-like_truncation_example") + trunc_inputs = TOML.parsefile(joinpath(trunc_dir, "gpec.toml")) + diiid(psihigh) = Equilibrium.setup_equilibrium( + Equilibrium.EquilibriumConfig(merge(trunc_inputs["Equilibrium"], Dict("psihigh" => psihigh)), trunc_dir), nothing) + long, short = diiid(0.9999), diiid(0.995) + + # One lower X-point, 1e-4 outside this file's boundary; the law reproduces the traced shear. + law = separatrix_q_law(short) + @test law !== nothing + @test law.psi_x[1] ≈ 1.0001 atol = 2e-5 + @test law.A[1] ≈ 0.640 rtol = 1e-2 + @test law.z_x[1] < 0 + @test law.shear_ratio ≈ 1.0 atol = 0.05 + @test separatrix_q_law(short; shear_ratio=1.001) === nothing + + # Anchored at psi = 0.995, it must place q = 7 and q = 8 where the grid to 0.9999 traces them. + for q_target in (7.0, 8.0) + psi_traced = find_zero(p -> long.profiles.q_spline(p) - q_target, (0.99, 0.9999), Bisection()) + gap_traced = law.psi_x[1] - psi_traced + @test law.psi_x[1] - separatrix_q_psi(law, q_target) ≈ gap_traced rtol = 0.05 + end + + # The scan's overlap point must therefore not depend on where the grid stops, and no + # rational may fall between the in-grid search and the law. + oi = GeneralizedPerturbedEquilibrium._overlap_slayer_inputs(trunc_inputs, trunc_dir) + one_n() = (i = ForceFreeStatesInternal(); i.nlow = 1; i.nhigh = 1; i) + scan_long = GeneralizedPerturbedEquilibrium._layer_overlap_scan(oi, one_n(), long) + scan_short = GeneralizedPerturbedEquilibrium._layer_overlap_scan(oi, one_n(), short) + @test scan_short.psihigh ≈ scan_long.psihigh atol = 5e-5 + @test scan_short.first_overlap == scan_long.first_overlap + @test all(diff(scan_short.m) .== 1) && all(diff(scan_long.m) .== 1) + @test any(scan_short.extrapolated) && any(occursin("separatrix law", note) for note in scan_short.caveats) + # psihigh stays the upper limit: at most one rational is placed past the equilibrium boundary, + # and the recommendation is never beyond it. + for scan in (scan_long, scan_short) + @test count(note -> occursin("skipped", note), scan.notes) <= 1 + @test scan.psihigh < 1.0 + end + + # No 2D flux, no X-point: an analytic equilibrium is not extrapolated. + sol_dir = joinpath(dirname(@__DIR__), "examples", "Solovev_ideal_example") + sol_inputs = TOML.parsefile(joinpath(sol_dir, "gpec.toml")) + solovev = Equilibrium.setup_equilibrium( + Equilibrium.EquilibriumConfig(sol_inputs["Equilibrium"], sol_dir), Equilibrium.SolovevConfig(sol_inputs["SOL_INPUT"])) + @test separatrix_q_law(solovev) === nothing +end