Repository navigation
Tearing.ModifiedRutherford - FEATURE! - Add a configurable steady modified Rutherford equation with seed thresholds, saturated widths and coil-forcing critical fields - #665
Open
logan-nc wants to merge 31 commits into
Conversation
…ford terms and F(w) output Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…diagnostics and format the submodule Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…e ModifiedRutherford sub-table Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
… linear gate and Callen preset The modified Rutherford solve now returns the seed threshold w_m and saturated width w_sat per surface (grid brackets refined by Brent), their minimum w_m_star, and validity flags (fully stabilized, unstable to any seed, resistive interchange, linear tearing, below w_ib, large Delta' w, complex Delta'). A [Tearing] layer solve in the same deck gates the roots and sets the lower end of the search from its layer width (API keyword `tearing`). The island frequency defaults to omega_*e (1 + eta_e/4) (omega_model), and delta_prime_factor, D_nc_factor and rho_theta_factor support scans. [Tearing.Callen] / CallenProblem(pe) adds Callen's reduced equation for RMP-driven islands on a coil-driven perturbed equilibrium: critical vacuum half-width, critical resonant field and the cubic roots, with the vacuum resonant field from the coil harmonic on each rational surface. The roots solve the cubic exactly, twice the Fortran 1.6.0 values. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
Callen 2016 Eq. 17 poses the steady cubic in the full island width with w_vac = 2h_v, so its exact roots are full widths. hw_min and hw_sat are now half of them, matching Fortran GPEC 1.6.0 (roots agree to 1.2e-3 on the DIII-D-like example at 10x coil current). Adds the Callen Fig. 8 roots as a unit test and a deck and API parity case for [Tearing.Callen]. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
… unstable surfaces, use true Delta' in Callen - Polarization rho_theta, w_ib and the layer width delta_r map metres to psi_N with the radial_label dr/dpsi used by SLAYER and LayerOverlap, replacing 2 psi_s / r_s. - Linearly unstable surfaces (gamma > 0) keep w_sat; w_m is NaN there. Only interchange instability or D_I >= 0 excludes a surface. - Callen preset takes Delta_c from the computed Delta' referred to rho_tor by default (main adds a Riccati solve); delta_prime_source = "cylindrical" restores Callen's -2m/rho. - The MRE accepts delta_prime_source = "cylindrical" (-2m referred to psi_N), so it runs from a Forward solve without Delta'. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…Rutherford equation saturation = true replaces r_s Delta' by r_s Delta'(0) - C_sat W/r_s (Fitzpatrick 2023 Eq. 9.15, after Hastie, Militello & Porcelli 2005), C_sat = 0.8 alpha_s^2 - 0.27 beta_s - 0.09 alpha_s, with alpha_s, beta_s from J_z = (2 - s)/q in the radial_label radius. beta_s takes r^2 (dimensionless), which reproduces Fitzpatrick's Fig. 9.2. The correction is referred to psi_N as the outer Delta' is and carried through the classical term's pressure factor. Off by default; the unit test checks the coefficients against the Wesson cylinder and the root against Eq. 9.21. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…franov with a windowed fit The cylindrical identity J_z = (2 - s)/q gives J_z <= 0 at q = 2 on shaped equilibria (the radial_label shear exceeds 2 inside the current channel), and its central differences need q''', which is piecewise constant on the spline, so C_sat did not converge. J_z is now (R0^2/B0) <mu0 j_phi/R> from the Grad-Shafranov flux-surface identity (current_density_profile), J_z' and J_z'' come from a least-squares quadratic over psi_N +- saturation_window (0.02, replacing saturation_dpsi), and saturation_unconverged flags surfaces whose C_sat moves more than 20% at half or double the window or whose window leaves the grid. On the TJ equilibrium J_z(0) = 2/q0 to 0.5% and C_sat at q = 2 matches the cylinder to 5%. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…rm matched to the SLAYER critical Delta Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…earingResult is saved after them The TearingResult writer now clears only its own part of Tearing/, and record_inputs! keeps recorded sub-tables a result does not state, so the saved file no longer depends on the order of the writes. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
… factor and fix review findings The Lütjens term was in Δ′_mn units while the other MRE terms are in ψ_N⁻¹; it now passes through classical_term as the saturation term does, so classical plus curvature change sign where SLAYER's gate does at every D_I (unchanged at D_I = −¼). large_delta_prime_w uses the dimensionless Δ′(w) w. Also: Callen single-n parse guard, the cubic at a = 0, the J_z sign probe, stale docstrings, and tests that reuse main's Riccati solve and drop assertions restating the code. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…ratio by their real keys Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…in's return Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…convention for D_I ≠ −¼ Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…d test header Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…ng root above the island start and flag unsaturated, off-grid and gate-disagreeing surfaces An island seeded at w_m grows to the first root above it, not the largest; where F(w_lo) > 0 it starts from w_lo. New flags: unsaturated (F > 0 at w_max), outside_plasma (psi +- w_sat/2 leaves the psi grid) and gate_disagrees (the sign of F(w_lo) against the linear growth rate). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…Callen module with output under Tearing/Callen/PerSurface Callen shares no geometry or terms with the modified Rutherford equation. Its result names the surface rational_psi and rational_q, as TearingResult does, and its datasets move to Tearing/Callen/PerSurface/. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…or every layer consumer layer_plasma returns the profiles, Z_eff, the main-ion mass and charge, the mass per electron, the resistivity closure and the chi splines once; the tearing solve, the layer-overlap scan, Callen, layer_eta and layer_rho use it in place of the SLAYER private helpers. Values passed to the layers are unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…d the result by keyword and drive the writer from its annotations The surface loop becomes _surface_inputs, _mre_terms, _gate and _surface_roots, and ModifiedRutherfordResult is built by keyword. The surface identity is rational_psi and rational_q, as in Tearing/PerSurface/; the writer emits every annotated field. The plasma comes from layer_plasma, so the ion mass is the one the layers use. Both Tearing subgroups are listed in hdf5-conventions.md. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…e the MRE harness case a variant of diiid_slayer_n1_toroidal Info/Runtimes gains callen_delta_prime. The Callen case reads Tearing/Callen/PerSurface/ and comments its override. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…h the seed threshold that reaches it When F crosses zero more than once, w_sat is the first of several saturated solutions. The result now keeps all of them (w_sat_all, branch-padded) with the rising crossing below each (w_seed_all), the seed above which the island grows past w_max (w_seed_unsaturated), the count n_sat and a multiple_saturated flag, and writes them under PerSurface. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…he vacuum island as wide as each seed threshold b_m, b_seed_all and b_seed_unsaturated are the resonant area-weighted fields [T] whose vacuum islands are as wide as w_m, each saturated branch's seed and the unsaturated seed, B = pi |n q'| chi1 w^2/(8A), with b_m_star their minimum. The new public vacuum_island_field inverts the perturbed equilibrium's island width; the API test checks it against that width and field at every surface. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…quation, with its critical fields and a Callen preset - Tearing.ResonantForcing: the resonant forcing of a coil-driven perturbed equilibrium (vacuum and total resonant fields, Callen's Delta_RMP, vacuum island width, toroidal-flux radius) shared by Callen and the MRE, and forcing_thresholds for the critical fields. - ModifiedRutherfordProblem takes pe; the forcing term Delta_RMP (drho/dpsi_N)(w_v/w)^2 enters F (forcing = false switches it off), so roots and branches include the applied field. b_crit (a forced window opens below the unforced seed threshold) and b_trigger (it reaches that threshold) follow from the unforced F. - delta_prime_model and pol_model gain :callen; CALLEN_PRESET (deck preset = "callen") reduces F to Callen's steady state, whose roots and b_crit equal CallenProblem's. - main shares one Riccati Delta' solve between the MRE and Callen when the force-free-states solve has none. - HDF5 output, docs, tests and the diiid_n1_callen harness case extended. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…allen preset replaces it The MRE's preset = "callen" reproduces Callen's cubic roots and critical field on the same resonant forcing, so the closed-form CallenProblem, its [Tearing.Callen] deck section and its Tearing/Callen/PerSurface output are removed; a deck with [Tearing.Callen] now stops at parse time. The closed form stays as the test reference for the preset, evaluated from resonant_forcing and Callen's polarization width, and the Callen 2016 Fig. 8 check now runs through the MRE root finder. main runs its Riccati Δ′ solve for the MRE alone, and the diiid_n1_callen harness case tracks the preset's quantities. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…fraction trapped_fraction used h = B_min/B_max and the inverted ratio <B>^2/<B^2>, which overestimates the effective trapped fraction by 25-40 % (0.7745 vs exact 0.6129 at epsilon = 0.2). It now returns the Lin-Liu & Miller 1995 upper bound f_tu with h = <B>/B_max and <B^2>/<B>^2, within 2 % of the exact value (above it, as an upper bound) for epsilon between 0.05 and 0.4. B_min is kept in the signature but is no longer used. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…Omega_E, normalize J_z by the vacuum field and flag collisional surfaces - omega_model = "exb" becomes "locked": an island at rest in the laboratory has omega = -n Omega_E in the E x B frame (the old option had the opposite sign). - current_density_profile normalizes by bt0 at rmean, the radius it is defined at, and stops with a clear error when the global parameters are missing. - New nu_star_e (Sauter 1999 Eq. 18b) per surface and a `collisional` flag (nu*_e > 1), since D_nc is the banana limit. - Docstring labels the Lutjens curvature away from D_I = -1/4 as this code's convention, and separates its W_d from the bootstrap w_d. - Tests: deck and API parity now covers every problem keyword (two variants), the saturation alpha_s, beta_s check uses closed forms in place of stored numbers, and the curvature test checks its large-w limit. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
…r's own keys, add Z_i to the poloidal Larmor radius and take B_pol² directly - The TearingResult writer clears the keys it writes under Tearing/ (an allow-list) instead of everything but a list of foreign groups, so any result can share Tearing/. - rho_theta = m_i v_Ti / (Z_i e <B_pol>); unchanged for hydrogenic ions. - B_pol² from the poloidal flux gradient instead of B² − (F/R)²; agrees to 2e-6 inside the plasma and is more precise near the edge. - The n_i and chi_e splines are built once per solve, not once per surface lookup. - Callen's polarization width is computed only for pol_model = "callen". Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae
logan-nc
marked this pull request as ready for review
November 14, 2026 02:47
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Adds a steady modified Rutherford equation (MRE) to
TearingasModifiedRutherfordProblem: every term switchable, seed thresholds and saturated widths per rational surface, and the coil forcing with its critical fields. Callen's reduced steady state becomes a preset of it, andCallenProblemis retired. The branch also carries a trapped-fraction fix inUtilities, previously its own PR.Closes #630. Part of #629: the critical resonant field and the vacuum-island widths are here; saturated-island pressure flattening is not.
What it computes
For each rational surface of a single-n force-free-states (FFS) solve, F(w) is the sum of the switched terms below, with w the full island width in ψ_N:
From F(w) it reports:
The deck section is
[Tearing.ModifiedRutherford](preset = "callen"for Callen's form); output goes toTearing/ModifiedRutherford/in gpec.h5. With[Tearing]present, the linear layer solve is the MRE's gate. Docs:docs/src/inner_layer.md.Also in this branch
Utilities.NeoclassicalResistivity.trapped_fractionevaluated a different expression from the Lin-Liu & Miller (1995) bound it cites. It ran 15–45 % above the exact effective trapped fraction and now evaluates Eq. (4) as written, to within 3 %. Its one caller isResistEval, which sets the F_33 neoclassical resistivity, so η falls and S rises on every surface. The recommended weighted average (Eq. 18) needs a further surface average and is left for later.layer_plasma(InnerLayer, new API): one view of the kinetic plasma for every layer consumer.Tearing/HDF5 output. ATearingResultwrite now clears only the keys it owns, so MRE output written before it is kept.Release note
runtests_equil.jl:312(TokaMaker i-file), which fails identically on develop.CallenProblemis removed; useModifiedRutherfordProblem(ffs; pe, CALLEN_PRESET...)or[Tearing.ModifiedRutherford] preset = "callen".New
[Tearing.ModifiedRutherford]section: a configurable steady modified Rutherford equation giving seed thresholds, saturated widths and coil-forcing critical fields per rational surface. The neoclassical trapped fraction behind the Sauter and Redl resistivity models was 15–45 % too high and now follows Lin-Liu & Miller Eq. (4), so η is lower and S higher on every surface.🤖 Generated with Claude Code
https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae