EQUIL - BUG FIX - Fix kappa calculation in separatrix finder#208
Merged
Conversation
…paratrix finder The zsep loop in equilibrium_separatrix_find! used unbounded Newton iteration to find dZ/dθ = 0 extrema. Newton diverged when d²Z/dθ² passed through zero (inflection point), sending θ to values like ±1e9. Due to sin/cos periodicity the diverged θ still produced plausible Z values, silently returning wrong extrema depending on floating-point luck. This caused kappa ≈ 0 or kappa < 0 for some equilibria. Fix: use bracketed Brent method on dZ/dθ with θ ∈ [0.5,1.0] (bottom) and θ ∈ [0.0,0.5] (top). This eliminates d²Z/dθ² entirely, is 4x faster, and guarantees convergence to the correct extremum. Adds regression tests for kappa with Solovev equilibria (e=1.6, e=1.0).
…force tests Convert the rsep loop (outboard/inboard midplane finding) from unbounded Newton to bracketed Brent for consistency with the zsep fix. Remove unused mtheta and vector variables that were only needed for Newton initial guesses. Add rsep value-level checks (amean ≈ 0.33, rmean ≈ 1.0) and circular cross-section rext symmetry test (rext[1] ≈ rext[2]).
logan-nc
approved these changes
Apr 2, 2026
Collaborator
|
Thanks!! |
logan-nc
added a commit
that referenced
this pull request
Apr 2, 2026
Incorporates: - PR #207: FastInterpolations v0.4 API (Series() wrappers, ZeroCurvBC, remove explicit defaults) - PR #208: Separatrix finder Brent method fix - PR #198: Analysis module improvements (plot_eigenvalues, plot_delta_prime, plot_ffs_summary) - AdaptiveArrayPools v0.3.5, FastInterpolations v0.4.7 Conflict resolved in src/Analysis/ForceFreeStates.jl: retained our plot_edge_stability_scan alongside develop's new plot_eigenvalues, plot_delta_prime, and plot_ffs_summary functions. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
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
2 participants
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.
Summary
equilibrium_separatrix_find!used unbounded Newton to find separatrix geometry. Without bracket constraints, Newton frequently diverged — θ escaped to values like ±1e2 or larger. Due tosin/cosperiodicity, the diverged θ still returned plausible values, but which extremum it landed on depended on spline coefficients and grid resolution. When both sides landed on the same extremum,kappa ≈ 0orkappa < 0.Fix
Replace all Newton calls in
equilibrium_separatrix_find!with bracketed Brent:dZ/dθ = 0, θ ∈ [0.5, 1.0] (bottom) / [0.0, 0.5] (top). Removesz_deriv2entirely.θ + η(θ) - η₀ = 0, θ ∈ [-0.25, 0.25] (outboard) / [0.25, 0.75] (inboard). Removes derivative-based Newton step.Adds regression tests for kappa and separatrix geometry with Solovev equilibria (e=1.6, e=1.0).
Performance (zsep loop,
@btime)No
d²Z/dθ²computation and fewer iterations within the bracket.