Repository navigation
Large non-hermitian component of plasma response matrix when calculating crit #93
Description
Activity
- changed the title
[-]Non-Hermitian plasma response matrix when calculated crit[/-][+]Large non-hermitian component of plasma response matrix when calculating crit[/+]on Dec 10, 2025 Progress update: I have been doing some investigating on this and will be adding some updates in an upcoming PR. I think this is an inherent issue in the serial EL integrator, where condition number of the solution matrix is becoming absolutely massive (and which isn't fixed by the fixups). After iterating with Claude a bit, it seems like the Ricatti path completely fixes this issue - since this seems to be what we are trending toward (and DCON has used this implementation for decades without too much problem), I think its ok to flag the bad conditioning when it appears (in the same place as the current non-Hermitian error). Unfortunately, I wasn't able to track down a full fix for this since I think its quite a complex issue.
To provide a bit more context from some of Claude's outputs that might help in the future:
Why is Wp Hermitian?
Wₚ⁻¹ = U₁·U₂⁻¹is built from the Euler-Lagrange fundamental solution matrixu(blocks
U₁ = u[:,:,1],U₂ = u[:,:,2]). Glasser (2016) —docs/resources/2016-Glasser-...pdf,
Eqs. 26/28/55/57/63/66 — showsWₚmust be Hermitian because of the symplectic invariant
U†JU = Jpreserved by the Euler-Lagrange ODE. Any measured anti-Hermitian part is pure
numerical noise, not physics. Physically, Wp is Hermitian because delta W_p is real.Why isn't it actually Hermitian?
The non-Hermitian warning fires on the serial standard Euler-Lagrange path when forming the plasma response matrix Wₚ⁻¹ = U₁·U₂⁻¹. I reconstructed U₁, U₂ from the HDF5 output and found that U₂ is severely rank-deficient at essentially every step — cond(U₂) ≈ 1e64 (median), up to 1e76, with ~28 of 35 singular values below machine precision near rational surfaces.
How to get rid of it:
- use_riccati = true fixes it directly for ideal fixed-boundary runs: it resets U₂ = I by construction, dropping cond(U₂) from ~1e64 to ~1 and nonherm max
from 1.6 to 0.06, with identical stability results (Newcomb crossing count preserved). - use_parallel = true (production default) avoids it too — per-chunk identity ICs keep the conditioning bounded. The issue is specific to the serial standard path.
However, this does not fully get rid of the problem, as there are still large spikes near the rational surfaces that appear due to the singular crossing methods.
- Riccati converts a large, structural, whole-domain problem into a small, localized, near-surface residual. That's a real improvement (the warning stops firing across the bulk of the domain), but it is not a complete fix, and the plot is honest about that.
- The remaining 0.14 is not removable by re-gauging U (QR, Riccati, etc.) because it isn't a U₂-conditioning artifact — it's the resonant-layer approximation.
Things that do NOT help
ucrit/Gaussian reduction (rebalances U₁ magnitudes, not U₂ rank), tighter tolerance, or sing_order. A right-gauge QR re-orthonormalization is mathematically valid but breaks the crossing bookkeeping (the singular-surface logic depends on the Gaussian-reduction index arrays it would suppress), so it's not a clean drop-in.
Changingsingfac_mincan help in some specific cases, but makes the Delta prime calculation blow up often and is unclear.Why we can't just always use Riccati
it's a serial path that can't produce the FM propagators needed for the free-boundary/STRIDE Δ′ matrix, doesn't support kinetic modes, and needs populate_dense_xi = true to feed PerturbedEquilibrium.
Suggested action
treat as cosmetic. Either point users to use_riccati = true when they hit it on the serial path, or gate the warning on cond(U₂) < ~1e14 so it only fires when the response matrix is actually trustworthy.
- use_riccati = true fixes it directly for ideal fixed-boundary runs: it resets U₂ = I by construction, dropping cond(U₂) from ~1e64 to ~1 and nonherm max
I think the TL;DR here is I am still unsure if this is an actual issue in the code or not and if this should remain an open issue. Either way, its been around in DCON for a while, and I can't tell if there's an actual mathematical basis for enforcing it or if we're hiding numerical issues - for example, I found in the DCON paper, JK's thesis, and Nik's thesis this symmetrization of the Wp matrix:
Glasser 2016
Park 2009
Logan 2015
@jhalpern30 did you explore if there was any additional Gaussian-reduction-like trick that could regularly fix the U2 rank and then be undone like the Gaussian reductions are? This reads like "Gaussian reduction only fixed half the problem and needs an extension in scope".
This had come up in the AI response - I think there's some merit to 1) below, but the mathematical justification required for that requires more human intervention than just "do it Claude" so I stopped there since I had already spent more time than I expected diagnosing this (but ended up learning quite a bit along the way). Specifically, I was concerned because Glasser focuses hard on the symplectic/Hermitian properties in the 2016 paper, and even says the below in Section III: Singular Surfaces, seemingly contradicting the finding that the resonant crossings break the matrix structure
But good point, I will document the suggestions here for future reference
- Post-crossing symplectic re-projection (highest value). cross_ideal_singular_surf!'s column-zeroing + trapezoidal jump + asymptotic substitution is the one step in the whole integration that isn't a symplectic transformation — it's what breaks U†JU=J. A cheap correction (project U back onto the constraint U₁₂†U₂₂ = U₂₂†U₁₂ right after the jump) would stop the defect from being injected in the first place, rather than relying on hermitianpart! to clean it up much later at every crit evaluation. This directly explains the earlier finding that a coarse crossing contaminates the entire outer region downstream.
- Higher-order crossing integration. Replace the 2-point trapezoidal Euler jump across dpsi with a few sub-steps of the existing Vern9 integrator. This attacks the same local truncation error without needing a smaller dpsi (which is what actually hurts Δ′).
- Extended precision at the crossing, not just in the BVP. extended_precision_bvp already promotes the parallel-path Δ′ solve to Complex{Double64}; the serial path's sing_get_ca/sing_get_ua asymptotic extraction in cross_ideal_singular_surf! doesn't currently get the same treatment. Extending it there would buy back digits without moving singfac_min at all.
I think in a follow-up session when I have the time/motivation I can look into this. I think the issue can be left open at least until then, since I do think it is a legitimate bug or at least there's a better treatment.
Reacted by Nikolas LoganAgreed. These look like reasonably promising things to try but it can be back burner for now
As discussed in #91, when calculating crit (the smallest eigenvalue of Wp), the original Fortran code forces Hermitian symmetry of the matrix. According to the DCON paper, this matrix should already be Hermitian by construction, so this was confusing. As part of this PR, we added in a warning if the non-Hermitian part of this matrix becomes too large, which I arbitrarily set at 1e-3; however, in some of the basic examples such as the DIIID ideal example, this warning is being raised. It seems like if this warning if being thrown on the simplest of examples, it likely should not be a warning and just be fixed or removed.
I think ultimately, we need to figure out if there being a large non-Hermitian component is a bug that we need to fix and are masking by forcing the Hermitian symmetry or if it not affecting the results and the warning should just be removed.