Vacuum - PERF - Serve conjugate toroidal residue classes from one factorization - #438
Draft
jhalpern30 wants to merge 2 commits into
Draft
Vacuum - PERF - Serve conjugate toroidal residue classes from one factorization#438jhalpern30 wants to merge 2 commits into
jhalpern30 wants to merge 2 commits into
Conversation
…operators via conjugate pairs to reduce the number of solves by around 2x. Again, this was all Claude and I need to look at this to clean it up
Contributor
|
This pull request is missing a reviewer. If you are not ready to name them, mark this pull request as a draft. |
…es, ~2.9x faster at nfp=2
The operator element type was chosen once per call rather than per mode family: without
stellarator symmetry it was ComplexF64 for every family whenever nfp > 1. That included the
self-conjugate families k=0 and, for even nfp, k=nfp/2, whose phases exp(-2πikd/nfp) = ±1 are
real and whose reduced operators are therefore real. The cause was the shared flat buffer, which
carries a single element type, so one general family forced complex on all of them.
Both element types are now carved from one Float64 backing store, the complex families through a
strided reinterpret (still a StridedArray, so LAPACK applies). Peak allocation is therefore
byte-for-byte unchanged when any family is complex, and halved when every family is real, which
is the case for nfp <= 2 and for any stellarator-symmetric run.
Selecting the type alone would have been a net 2.6x regression: BLAS has no mixed real/complex
triangular solve, and Julia's generic fallback is ~30x slower than a matched one, costing more
than the real factorization saves. A real family now carries its complex right-hand side through
both the projection and the solve as a real [Re Im] pair, which is exact because a real operator
propagates the two parts independently. That path is 30x faster than the generic fallback and
1.9x faster than the matched complex solve, so this also removes the same latent penalty from the
existing nfp == 1 and stellarator-symmetric paths, which were already building real operators.
Measured, 8 threads, no wall, non-stellarator-symmetric boundary, min of 3:
nfp=2 N=9216 Np=4608 2.891 s -> 1.003 s 2.9x
nfp=3 N=9216 Np=3072 1.066 s -> 0.707 s 1.5x
nfp=4 N=9216 Np=2304 0.800 s -> 0.526 s 1.5x
nfp=3 N=9216 stell-sym 0.476 s -> 0.286 s 1.7x
The factorization itself gains ~3.5x; end-to-end is diluted by assembly, so the speedup grows
with resolution (nfp=2 goes 1.2x -> 2.2x -> 2.9x over N = 2048, 4608, 9216). nfp=2 gains most
because every family there is self-conjugate; it was previously the case where the field-period
reduction bought no factorization speedup at all. Peak RSS at Np=6144 drops 2.57 GB -> 2.23 GB.
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
1 participant
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.
Release note
Toroidal residue classes
kand-know share one operator assembly and one LU factorization, since the field-period blocks are real and soD̂₋ₖ = conj(D̂ₖ). A 3D stellarator run whose mode set spans conjugate classes gets 1.55-1.95x off the vacuum solve for no change in results.Regression report
regress --cases diiid_n1,solovev_n1 --refs feature/vacuum-3d-inductance,performance/vacuum-3d-conjugate-classes@185acd08diiid_n1solovev_n1Both cases are
nfp = 1, which has exactly one self-conjugate class, so no pairing can occur in either — they confirm nothing regressed, not that the reduction works.Notes for reviewers
Base branch is
feature/vacuum-3d-inductance, notdevelop— this stacks on that PR and should merge into it.Review package (diagrams, derivations, benchmarks): https://claude.ai/code/artifact/17fac702-393d-4993-9b49-e290234d3ecc
It covers this change and the branch below it as one story.
The identity
The field-period blocks
D_d,S_dare real for any boundary — that is already why the code picksFloat64whennfp == 1. SoD̂₋ₖ = conj(D̂ₖ)andŜ₋ₖ = conj(Ŝₖ), and withF := conj(E')for the partner class's Fourier rows:So the partner needs no new code path: run the existing per-class body with
conj(E')as its mode basis against the representative's factorization, then conjugate the output block. That expression is invariant under a simultaneous unitary change of basis, which is why the same rule works verbatim on the stellarator-symmetric path — the class-kbasis is reused and the real transformed operator is the same matrix for both classes.I_vis accumulated pre-conjugation and the loop already applies a trailingconj!, so for a partner the two cancel and that call is skipped. That one-line asymmetry is the likeliest place for a sign error —Vacuum.jl, thepartner || conj!(Iv_block)line — and is what the A/B matrix exists to catch.Verification
nfp1-5, wall/nowall,use_symmetryon/off,compute_Ivon/off)nfp=1,nfp=2, single-class mode family (no pair possible)test/runtests_vacuum.jlexamples/Solovev_ideal_example_3Dvs base branch, 135 numeric datasetsAgreement is roundoff rather than bitwise on purpose: the unpaired path builds
cis(-2*pi*(nfp-k)*d/nfp)while the paired path conjugatescis(-2*pi*k*d/nfp), and those differ in the last bits.Measured
nfp=5, 64x32 per period, conformal wall, 6 threads, back-to-back in one process:n=[2,3]n=1:4n=0:4The two mechanisms are multiplicative and do not interfere. Memory is neutral — peak is still one class's operator.
Worth pushing on
nfp = 1, which has exactly one, self-conjugate, class. The gain is entirely prospective.nfp) or 2 (evennfp) unpairable self-conjugate classes, so the realistic band is 1.3-1.6x; 1.79x needs a mode set that avoidsk=0andk=nfp/2.use_symmetrynoruse_conjugate_pairingis reachable from a full GPEC run. All three callers outside the Vacuum module (Free.jl:105,Free.jl:240,ResponseMatrices.jl:304) omit them. I kept them as kwargs rather thanVacuumInputfields becauseVacuumInputstates the physics problem and these do not change the answer — but if you want an escape hatch fromgpec.toml,ForceFreeStatesControlis the natural home. Not done here; it pairs naturally with plumbingnfpthrough.Commit message
The single commit subject is informal and says the work is unreviewed. Worth rewording before merge.