feat(lpse2d): 2-D hybrid particle evolution, exact-rotation light coupling, TPD pump depletion - #368
feat(lpse2d): 2-D hybrid particle evolution, exact-rotation light coupling, TPD pump depletion#368physicistphil wants to merge 2 commits into
Conversation
Lift the `ny == 1` restriction on `terms.hpe`. The Follett-2017 test-particle
Landau-damping feedback now runs in a genuinely 2-D box, with quasi-1D falling
out as the `n_angles == 1` special case of the same code path.
The physics: Eq. 4 integrates d<F>/dv over the resonance surface omega = k.v,
which is a point in 1-D but a line in (vx, vy). That line integral is exactly
the 1-D formula applied to the projection of f onto the mode's own direction,
P_khat(v) = int d^2v' f(v') delta(v'.khat - v)
gamma_L(k) = -(pi/2) (wp0^3/|k|^2) dP_khat/dv |_{omega_res/|k|}
so the extraction bins particles along `n_angles` directions spanning [0, pi)
and interpolates between the two bracketing directions per mode. Opposite
directions are mirror images, which is precisely what the existing sgn(kx)
encoded -- the 1-D sign convention was the d=1 shadow of this bookkeeping.
Implementation notes:
- angle/velocity gather indices are static (v_phi comes from the k grid), so the
per-step extraction is 4 gathers per mode from a flat (n_angles+1, nv) table,
with no (n_modes, nv) intermediate
- all angles are binned in one fused scatter, which keeps the cost flat in
n_angles rather than linear
- the tail cut is on the speed in 2-D; the projection of that speed-truncated
Maxwellian equals the true 1-D Maxwellian exactly above the cutoff, i.e.
everywhere the resonance reads it. tail_fraction() is now dimension-aware and
shared with the fhot diagnostics.
- thermalizing x wall gains a cosine-weighted emission direction and a tabulated
s^2 exp(...) speed law in 2-D; the 1-D closed form is untouched
- `y_thermal_frac` thermalizes a fraction of transverse *crossings*, counted in
the substep loop where the wrap happens
- state is (Np, ndim); epw_hist is (n_angles, nv) and is saved/plotted per angle
Also fixes, all latent while ndim == 1:
- fhot/hpe_mean_energy used a per-component u instead of |u|
- the smoke tests asserted every series value finite, contradicting the
documented NaN sentinel that hpe_gamma_ratio_kpeak emits before the EPW has
band energy (this failed on lpse2d/srs before this branch)
- the same tests asserted a single-mode damping ratio in one sample; kpeak is
shot-noise-limited per mode (frac < 0.5 is 0.039 at 20k particles, 0.016 at
200k) and the EPW noise seed is redrawn per run, so they now take a median
over the tail of the run
Found while validating the 2-D smoke run, and mitigated: the rate is read from
df/dv at one velocity bin per mode, so with finite particles some modes draw a
slope steep enough that the gamma >= 0 clamp sends them to exactly zero -- and a
mode pinned at zero is *undamped*, so it grows relative to the band instead of
the error averaging out. The clamp selects for its own errors. 2-D meets this far
more often (6644 band modes vs 668 in the smoke configs). New
`terms.hpe.hist_smooth` (binomial filter over the histogram before the slope,
default 2 in 2-D and 0 in quasi-1D so existing 1-D results stay reproducible) is
the mitigation, and is bias-free because the per-k calibration is derived through
the same operator. Fraction of band modes clamped to zero at hist_smooth 0/1/2/4:
0.061/0.035/0.019/0.004 at 20k particles, 0.0075/0.0021/0/0 at 200k, with the
band-mean ratio staying at 1.02-1.06.
Validation: tests/test_lpse2d/test_hpe_2d.py (11 tests). The load-bearing ones are
that an isotropic tail gives a rate matching the analytic Landau rate with no
angular bias, and that the ky = 0 column of a 2-D box reproduces the quasi-1D
result -- i.e. this is a strict generalization. Also covered: the angle-
interpolated dP/dv against a directly binned line integral for off-axis modes,
2-D loading/free-streaming/transverse wrap, the Ey gather branch, transverse
thermalization, and the absorbing-x wall in both geometries.
Cost (CPU, 400x120, Np = 5e5): the 2-D add-on is ~8x the 1-D add-on and flat in
n_angles (386/357/402 ms at 16/32/64). Scaled by the measured 3090 1-D add-on
that is roughly 2-4x a fluid run. n_angles is an accuracy knob, not the cost
knob, and 2-D needs no particle-count increase over 1-D.
Design record and measurements in docs/dev/lpse2d-hpe-plan.md Sec. 7; reference
docs in docs/source/solvers/lpse2d/config.md; production config in
configs/envelope-2d/srs-hpe-2d.yaml.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…t filter Three opt-in terms.light options for the coupled (pump-depletion) light solver, all default off so existing configs are unchanged: - coupling: rotation -- Strang-split each light sub-step as [exact E0<->E1 coupling over dt_l/2] [staggered propagation with the coupling off] [exact coupling over dt_l/2]. With laplacian phi frozen the coupling matrix satisfies M^2 = -Omega^2 I, so exp(tau M) = cos(Omega tau) I + sin(Omega tau)/Omega M is applied exactly and conserves w1|E0|^2 + w0|E1|^2 pointwise. The MATLAB staggered real/imaginary update is a leapfrog only for i x (real operator); on the complex, spatially rotating laplacian phi it is forward Euler on a rotation and grows at ~0.04 Omega per unit time independent of dt_l, which is what drove the 2D SRS+TPD testbed runs non-finite once delta n/n ~ 0.3. - tpd_depletion: true -- add i e/(4 wp0 me) exp(i (w0 - 2 wp0) t) rho E_y to dE0_y/dt, the energy-conserving partner of the TPD potential source (the MATLAB original depletes the pump through SRS only, so the EPW energy could exceed the absorbed laser energy). - filter: <fraction of pi/dx> -- optional isotropic low-pass on E0, E1 once per EPW step (diagnostic). RamanLight.rhs gains a couple flag so the E1 RHS can drop the exchange term when the rotation handles it. Tests: explicit-scheme growth for a complex laplacian, rotation action conservation to 1e-12, small-angle agreement with the explicit scheme, and dW_pump == -dW_epw for the TPD pair. Evidence: srs-campaign sims/lpse/srs-2d-testbed/NOTES.md, "Blowup investigation" (MLflow experiments lpse_2d_testbed / lpse_2d_testbed_blowup). Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Comparison with the 2D HPE implementation currently on
|
| Area | Current main / PR#363 |
PR#368 / e7b0903 |
|---|---|---|
| Angular representation | Oriented directions over [0,2pi) |
Independent axes over [0,pi) plus sign/mirroring |
| 1D/2D structure | Separate pusher, boundary, and damping branches | Unified (Np, ndim) state and shared operator |
| Histogram | lax.map of one histogram per direction |
Matrix projection plus one fused scatter |
| Damping interpolation | Angle/velocity indices formed during extraction | Static four-point angle/velocity gather indices |
| Smoothing | None | hist_smooth, default 2 passes in 2D |
| Initial positions | Full flattened background_density[ix,iy] |
background_density[:,0] for x; uniform y |
| Transverse boundaries | Periodic or thermalizing y, including corner handling | Periodic y only, optionally rethermalizing a fraction of crossings |
| Tail cap | Properly truncated Rayleigh draw below 0.99c |
Unbounded draw followed by minimum(..., 0.99c) |
| Quasi-1D raw state | Historical (Np,), (nv,) shapes |
(Np,1), (1,nv) internally |
| Default angular resolution | 32 directions over 2pi = 16 independent axes |
64 independent axes over pi |
What should be ported from e7b0903
[0, pi)projection axes- signed resonance and mirrored angular closure
- dimension-generic resonance/damping operator
- precomputed angle/velocity gather indices
- fused histogram implementation, after GPU peak-memory measurement
- exact-bin expected projected Maxwellian
y_thermal_frac- axial-reduction, direct-projection, smoothing, and transverse-crossing tests
What should be retained from current main
- correct modular periodic field gathering
- full 2D density-weighted position loading
- cap-safe tail sampling and associated input checks
- independent x/y boundary handling and corner behavior
- existing state/artifact compatibility
- current TPD/SRS/IAW integration
- the synthetic directional test that flattens one projection and verifies that only aligned modes are undamped
Histogram performance and scaling
The fused histogram is a substantial implementation improvement. In an illustrative local CPU benchmark with 100,000 particles, I measured:
PR#368 fused histogram
8 axes: 3.9 ms
16 axes: 8.6 ms
32 axes: 26.5 ms
64 axes: 69.0 ms
current main, full-circle lax.map
16 directions: 279 ms
32 directions: 497 ms
64 directions: 798 ms
These are not production GPU measurements, but they show that the fused implementation is much faster on this backend.
The description that its cost is "flat in n_angles" should nevertheless be qualified. It forms and processes (Np, n_angles) projected samples and indices, so histogram work and temporary storage still scale with Np * n_angles. The total HPE step can look approximately flat because hundreds of particle substeps dominate. At Np=500,000, n_angles=64, the projected-value and index temporaries may be material on GPU; peak memory and runtime should be measured there, with chunking considered if necessary.
Also, the meaning of n_angles changes. Main's default 32 directions over 2pi represents only 16 independent axes, whereas this PR's default 64 over pi represents 64 independent axes—four times the unique angular resolution. A directly comparable setting would be about 16 axes here. The default should be chosen by nonlinear angular convergence rather than silently inherited.
Histogram smoothing
The clamp-selection problem identified here is real: noisy positive slopes are clamped to zero damping, and accidentally undamped modes preferentially grow. hist_smooth is a useful mitigation.
However, per-k calibration makes smoothing unbiased only for the calibration distribution:
smooth_operator[f_Maxwell] -> analytic gamma_L
It does not make the smoothed operator identical to the unsmoothed operator for an evolved plateau or sharp tail feature. Smoothing changes nonlinear velocity resolution. Please document it as estimator regularization and require convergence over hist_smooth, nv, n_angles, and particle count. Logging a non-applied unsmoothed diagnostic alongside the smoothed applied rate would also be useful.
The statement that no particle-count increase over 1D is needed should also be softened. Every projection receives all particles, so particles are not divided among angles, but the many-mode extreme-value problem remains. The 2D smoke test itself raises n_particles from 20,000 to 200,000 because thousands of resonant modes otherwise produce clamp-selected zero-damping outliers.
Validation performed
I ran the non-slow HPE tests on both versions:
current main: 10 passed, 3 deselected
PR#368: 17 passed, 4 deselected
The sign/mirror bookkeeping also passed a synthetic anisotropic check: flattening the x projection undamped both +x and -x modes while leaving +/-y modes at the analytic rate.
The new test suite is stronger in several areas, especially axial reduction and direct off-axis projection. It should be combined with main's directional-flattening and full boundary tests. A quantitative nonlinear oblique-wave trapping test is still needed; the isotropic Maxwellian tests primarily validate geometry and the calibrated linear limit.
Proposed resolution
Use current main as the integration base and port the mathematical/operator improvements from e7b0903, rather than merging this branch wholesale. This PR is currently merge-conflicted with main and also bundles a separate light-coupling/TPD-depletion commit that overlaps work already merged through PR#363; the HPE and light reconciliations should be split.
Bottom line: PR#368 has the better fundamental mathematics; current main has the safer implementation. The target should be PR#368's half-plane/unified resonance operator inside main's density loading, periodic gather, boundary, compatibility, and integration framework.
…lter Two opt-in terms.light options for the coupled (pump-depletion) light solver, both default off so existing configs are unchanged: - coupling: rotation -- Strang-split each light sub-step as [exact E0<->E1 SRS exchange over dt_l/2] [staggered propagation with the exchange off] [exact exchange over dt_l/2]. With laplacian phi frozen the exchange matrix satisfies M^2 = -Omega^2 I, so exp(tau M) = cos(Omega tau) I + sin(Omega tau)/Omega M is applied exactly and conserves w1|E0|^2 + w0|E1|^2 pointwise. The MATLAB staggered real/imaginary update is a leapfrog only for i x (real operator); on the complex, spatially rotating laplacian phi it is forward Euler on a rotation and grows at ~0.04 Omega per unit time independent of dt_l (Omega = e|laplacian phi|/(4 me sqrt(w0 w1))), which is what drove 2D SRS+TPD pump-depletion runs non-finite once delta n/n ~ 0.3. The TPD pump term and the IAW detuning stay in the staggered RHS under both settings; with SRS off there is no exchange and the settings coincide. - filter: <fraction of pi/dx> -- optional isotropic low-pass on E0, E1 once per EPW step (diagnostic). RamanLight.rhs and CoupledLight.pump_rhs gain a `couple` flag so the RHS can drop the exchange when the rotation handles it. terms.light is now validated through LightModel in get_solver_quantities (as terms.iaw/hpe already are), and the two options are refused without pump_depletion. Tests (tests/test_lpse2d/test_light_coupling.py): explicit-scheme growth 1 + sin^2(arg L)(Omega dt_l)^2/2 for a complex laplacian, rotation action conservation to 1e-12, small-angle agreement with the explicit kick, rotation == explicit with phi = 0 over a full EPW step, action conservation of the rotation halves over a full sub-step loop, E1 stays zero with SRS off, filter mask semantics, and config validation. Evidence: srs-campaign sims/lpse/srs-2d-testbed/NOTES.md, "Blowup investigation" (MLflow experiments lpse_2d_testbed / lpse_2d_testbed_blowup). Light-only successor to #368. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01U2qwF9isiJGwvF5y6mbeGC
1 participant
Two commits on
lpse2d/srssince #320.1. 2-D hybrid particle evolution (e7b0903)
Lifts the
ny == 1restriction onterms.hpe. The Follett-2017 test-particle Landau-damping feedback runs in a genuinely 2-D box; quasi-1D is then_angles == 1special case of the same code path.n_anglesdirections in [0, π) and each mode interpolates between its two bracketing directions.n_angles; no particle-count increase over 1-D is needed.y_thermal_fracthermalizes a fraction of transverse crossings.terms.hpe.hist_smooth(binomial filter on the histogram before the slope; default 2 in 2-D, 0 in quasi-1D). Motivation: the γ ≥ 0 clamp pins shot-noise-steep modes to exactly zero, and a zero-damped mode grows relative to the band. Fraction of band modes clamped at hist_smooth 0/1/2/4: 0.061/0.035/0.019/0.004 at 20k particles.fhot/hpe_mean_energyused a per-component u instead of |u|; HPE smoke tests now honour the documented NaN sentinel ofhpe_gamma_ratio_kpeakand take a median over the tail of the run.tests/test_lpse2d/test_hpe_2d.py(11 tests: isotropic tail matches the analytic Landau rate with no angular bias; the ky = 0 column of a 2-D box reproduces quasi-1D) plus the updatedtest_hpe.py. Design record indocs/dev/lpse2d-hpe-plan.md§7; configconfigs/envelope-2d/srs-hpe-2d.yaml.2. Exact-rotation light coupling, TPD pump depletion, light filter (abe1b9e)
Three opt-in
terms.lightoptions for the coupled (pump-depletion) light solver, all default off:coupling: rotation. The MATLAB staggered real/imaginary light update is a leapfrog only for operators of the form i × (real operator). The E0↔E1 coupling through the complex, spatially rotating ∇²φ is not one, so on that term the scheme is forward Euler on a rotation: per-sub-step amplification 1 + sin²(arg ∇²φ)(Ω dt_l)²/2 with Ω = e|∇²φ|/(4mₑ√(ω₀ω₁)), and with propagation in the loop a growth of about 0.04 Ω per unit time independent of dt_l. With ∇²φ frozen over the sub-steps the coupling matrix satisfies M² = −Ω²I, so exp(τM) = cos(Ωτ)I + sin(Ωτ)/Ω·M is applied exactly, Strang-split around the existing propagation update, and conserves ω₁|E0|² + ω₀|E1|² pointwise.tpd_depletion: true. Adds i e/(4ωₚ₀mₑ)·e^{i(ω₀−2ωₚ₀)t}·ρE_y to dE0_y/dt, the energy-conserving partner of the TPD potential source (the MATLAB original depletes the pump through SRS only, so the EPW energy could exceed the absorbed laser energy).filter: <fraction of π/dx>. Optional isotropic low-pass on E0, E1 once per EPW step (diagnostic).RamanLight.rhsgains acoupleflag so the E1 RHS can drop the exchange term when the rotation handles it. Tests intests/test_lpse2d/test_light_coupling.py: explicit-scheme growth for a complex Laplacian, rotation action conservation to 1e-12, small-angle agreement with the explicit scheme, dW_pump = −dW_epw for the TPD pair.Evidence from controlled runs (srs-campaign
sims/lpse/srs-2d-testbed/NOTES.md, "Blowup investigation"; MLflowlpse_2d_testbed/lpse_2d_testbed_blowup): the first 2-D SRS+TPD testbed runs went non-finite at 4–5 ps with the explicit scheme (halving dt_l only moved the NaN from 4.50 to 4.85 ps);coupling: rotationruns stay finite to 20 ps; with TPD on but no depletion the EPW held 4× the absorbed laser energy at 4 ps, andtpd_depletionholds that ratio at ≤ 0.53.Behaviour of existing configs
Unchanged: quasi-1D HPE runs use
hist_smooth: 0and then_angles == 1path; the three light options default off.Test status
tests/test_lpse2d/test_light_coupling.py: 4 passed (CPU, laptop).tests/test_lpse2dsuite (pytest tests/test_lpse2d -x, CPU jax 0.9.0.1, laptop): 60 passed, 1 skipped, 4 deselected in 5 min 35 s.Not in this PR
coupling: rotationandtpd_depletionshould become defaults. The rotation scheme still has narrow resonant windows above Ω·dt_l ≈ 0.009 in the periodic model problem; a unitary (spectral or Crank–Nicolson) light propagator would remove them.🤖 Generated with Claude Code