fix: Vlasov1d vth conventions#311
Merged
joglekara merged 12 commits intoJul 20, 2026
Merged
Conversation
The engine convention for _vlasov1d (and _vlasov2d, _pic1d) is the RMS/standard-deviation thermal speed: Maxwellian exp(-v^2/2) at T=1, L0 = lambda_De, code wavenumber k*lambda_De. normalization.py:91 was the sole outlier with v0 = sqrt(2 T0/m_e), which made every logged physical unit (v0, x0, c_light, box_length) wrong by sqrt(2), every dimensional string input (um box sizes, gradient scale lengths, laser k0) off by sqrt(2), and a 'normalizing_temperature: 2000eV' run physically a 4000 eV plasma. Dimensionless numeric-input dynamics were unaffected. Also fixes a sign error in the NRL Coulomb logarithm (log(n^0.5 / T^-1.25) == log(n^0.5 * T^+1.25)), which produced logLambda_ee = -11.8 and a negative logged nuee at 2000 eV / 1.5e21cc; the correct values are +7.22 and a positive rate. vth_norm() is deliberately untouched: its only callers are in vfp1d, which is self-consistently built on the sqrt(2T/m) convention. Adds tests/test_vlasov1d/test_units_boundary.py: dimensional input in, dimensional quantity out, checked against CODATA constants and the NRL formulary - the units-boundary test class whose absence let this bug survive every existing suite. See VLASOV1D_CONVENTIONS_AUDIT.md (F1, F2) and ADEPT_CONVENTIONS_AUDIT.md. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…in write_units _tf1d re-implements the debye normalization inline; it carried the same v0 = sqrt(2T/m) outlier (its engine, like _vlasov1d, runs in sqrt(T/m) units) and the same log(n^0.5 / T^-1.25) sign error. Diagnostics-only: the logged units were wrong; dynamics are unaffected (grid beta is written but never consumed). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…T0/mass - Dougherty.compute_vbar returned the raw first moment n*u while its callers (find_self_consistent_beta, the drag term) treat it as a mean velocity. The drag therefore centered on n*u and momentum was not conserved wherever n(x) != 1. Now returns sum(v f)/sum(f). - The Krook target Maxwellian was hard-coded to exp(-v^2/2) (T=1, m=1): any species with T0 != 1 was dragged toward T=1. It is now built with variance T0/mass from the species' bulk parameters, which are threaded through cfg.grid.species_params (new T0 entry). - New test drives the collision operator in isolation with a 50% density modulation and a finite drift: per-cell density, momentum, and energy conservation, and the relaxed mean velocity must equal u0 (the old operator drifts it toward n*u0). Audit findings F4, F5. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…rgy monitor - Field moments p and q were centered on the raw first moment n*u; they are now centered on the true mean velocity. The first moment is saved under the honest name 'j' and 'v' is now j/n (previously 'v' held the flux while the same integrand was called 'mean_j' in the scalars). - The field-moment '-flogf' computed +int f log|f| dv, the exact negative of the scalar of the same name; both now agree. - Adds mean_kinetic_energy / mean_field_energy / mean_total_energy scalars (electrostatic energy monitor, with the 1/2 factors). A conservation monitor of this kind would have caught the sqrt(2) convention bug far earlier. Audit findings F6, F7, F12.2. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
dx was xmax/nx and the interpolation period was xmax; both are now the box length (xmax - xmin). Latent for all shipped configs (xmin: 0.0) but wrong spacing, kx grid, and advection wrap for any xmin != 0. Audit finding F11. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
GridConfig.c_light (set only by wavepacket.yaml) was never read - beta/c_light are always derived from the normalization. The AmpereSolver class docstring claimed j = sum_s (q_s/m_s) int v f dv; the code correctly has no 1/m_s. Audit findings F12.1, F12.4. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The ey driver pair was a backscatter-matched triad under the old c_hat = 11.30 (w_pump - w_seed = 1.16 = old EPW frequency) with placeholder k0 = 1.0 on both drivers. Retuned to a matched triad under the corrected c_hat = 15.984: pump (k0 = +0.162949, w0 = 2.79), seed (k0 = -0.086080, w0 = 1.700942), EPW at k*lambda_De = 0.249. Matching condition and conventions documented in the config. Audit finding F3. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
- config.md gains a 'Normalization convention' section: v0 = sqrt(T0/m_e) (sigma convention), L0 = lambda_De, k in 1/lambda_De, Maxwellian exp(-v^2/2) at T=1, Bohm-Gross w^2 = 1 + 3k^2. - Species v0/T0 documented as code-units-only (they bypass normalize()); drift and thermal width now share the same unit. - FP/Krook 'baseline' documented as a rate in units of wp0 (no 2pi), with the newly logged nuee_norm as the reference scale, and the O(1) caveat between the Dougherty nu and the NRL nu_ee. - Super-Gaussian alpha documented (code comment + config.md): it fixes <v^4>/<v^2> = 3 T0/mass for all m; the variance equals T0/mass only at m=2 (x1.24 at m=3, x1.37 at m=4). Documentation only, per review. Audit findings F8, F9, F10 (nuee_norm logging itself landed with the species_params change in an earlier commit). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The fixtures locked in the sqrt(2)-wrong units (c_light 11.302, v0 2.652e7 m/s, x0 12.14 nm, negative nuee). Regenerated under the corrected normalization and the new config surface: c_light 15.984, v0 1.876e7 m/s, x0 8.584 nm, logLambda_ee +7.22, positive nuee, new nuee_norm, species_params T0, c_light knob removed. Co-Authored-By: Claude Fable 5 <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
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.
(Claude-generated, with some light edits)
Fix vlasov1d normalization convention (√2 in v0) and related conventions audit findings
Summary
Two independent conventions audits (
2026-07-16_VLASOV1D_CONVENTIONS_AUDIT.mdand2026-07-16_ADEPT_CONVENTIONS_AUDIT.md, now removed) found thatelectron_debye_normalizationinadept/normalization.pyused the most-probable-speed conventionv0 = sqrt(2*T0/m_e), while the entire_vlasov1dengine (initializer, collision operators, gold tests, all numeric configs) is built on the RMS conventionv0 = sqrt(T0/m_e). This branch fixes that root cause plus the other confirmed findings from the audits, scoped to_vlasov1dand_tf1d.The dimensionless dynamics were never wrong — the solver core never sees
v0, so numeric-input runs are byte-for-byte unchanged. What was wrong is the physical interpretation bolted on at the unit boundary:normalizing_temperature: 2000eVandT0: 1.0was physically a 4000 eV plasma.L0wassqrt(2)·λ_Deinstead ofλ_De, so every dimensional-string input (box sizes in µm, gradient scale lengths, laserk0) was off by √2, and every logged unit (v0,x0,c_light,box_length) was √2-wrong.logLambda_ee = -11.8(and a negative loggednuee) at 2000 eV / 1.5e21 cm⁻³; corrected to +7.22 and a positive rate.Changes
Core fix
adept/normalization.py:v0 = sqrt(T0/m_e)(RMS convention), sox0is the true Debye length; fixed the Coulomb-log sign (log(sqrt(n0) * T0_eV^-1.25)).vth_norm()is deliberately untouched — its only callers are invfp1d, which is self-consistently built onsqrt(2T/m)(sentinel:test_kappa_eh).adept/_tf1d/modules.py: same two fixes in tf1d's private inline copy ofwrite_units(diagnostics-only; tf1d dynamics don't consume these).Collision operators (audit F4, F5)
compute_vbarwas returningn·ū(missing1/n) — momentum was not conserved wherevern(x) ≠ 1. Nowsum(v·f)/sum(f).exp(-v²/2), dragging anyT0 ≠ 1species toward T=1. It now uses the speciesT0/mass(threaded throughcfg.grid.species_params). This path was latent (is_on: Falsein all shipped configs).Diagnostics (F6, F7)
"v"was actually the flux∫v f dv; it is now saved honestly as"j", and"v"holds the true mean velocityj/n. Pressure and heat-flux moments are now centered onū, notn·ū.-flogfentropy had the wrong sign (disagreed with the scalar version); fixed.nuee_norm(code-unit collision rate, the reference forbaseline),mean_kinetic_energy,mean_field_energy,mean_total_energy(electrostatic energy-conservation monitor).Grid (F11, latent)
dx = (xmax - xmin)/nxand the semi-Lagrangian interpolation period isxmax - xmin; both previously ignoredxmin. No shipped config usesxmin ≠ 0.Config / cleanup (F3, F12)
GridConfig.c_lightknob (always derived from the normalization; onlywavepacket.yamlset it, never read).srs-debug-small.yamlnumeric ey drivers to the correctedĉ = 15.984: pumpk0 = 0.162949, seedk0 = -0.086080,w0 = 1.700942, EPW atkλ_De = 0.249, with the backscatter matching condition documented inline. (Under the oldĉ = 11.30the triad was internally inconsistent between the ES and EM branches.)AmpereSolverdocstring (claimed a spurious1/m_sin the current).Docs & tests
docs/source/solvers/vlasov1d/config.mdnow states the normalization convention explicitly (v0 = √(T0/mₑ), L0 = λ_De, Maxwellianexp(-v²/2)at T=1, Bohm–Grossω² = 1 + 3k², ĉ ≈ 15.98 at 2 keV) and documents the unit traps: speciesv0/T0are code-units-only,baselineis a rate in units of ω_p0 (no 2π) withnuee_normas reference, and the super-Gaussian α pitfall (variance = T0/mass only at m=2; ×1.24 at m=3, ×1.37 at m=4).test_units_boundary.py(dimensional round-trip against CODATA + NRL — the class of test whose absence let the √2 bug survive) andtest_fp_momentum_conservation.py(isolated Dougherty operator with 50% density modulation + drift; per-cell n/momentum/energy conservation).c_light 11.302 → 15.984,x0 12.14 → 8.584 nm,logLambda_ee -11.78 → +7.22.Breaking changes / migration
2000eV" run was really 4000 eV. Anyone who back-computed physical parameters from old runs, or tunednormalizing_temperatureto hit a target, must re-interpret.srs.yamlwithxmax: 100um): these now map to the correct code values, so results will differ from pre-fix runs."v"as a flux must switch to"j"; the-flogfsign also flipped.grid.c_lightmust remove it (knob deleted; it was dead).srs-debug-small.yamlor other numeric EM drivers tuned againstĉ = 11.30need retuning toĉ = 15.984.nueeshould note the sign/magnitude correction.Out of scope
The audits' findings for other solvers are not in this branch:
_vlasov2dKrook target,_pic1d(inherits F10/F11),vfp1dinverse bremsstrahlung,lpse2d, and the orphanedvlasov1d2v.vth_norm()is intentionally unchanged to protectvfp1d.Test plan
Run locally on 2026-07-17 after the branch commits. Only pre-existing failures unrelated to this branch remained (3 tf1d resonance tests broken on main via
KeyError: 'ax', plus osiris deck-roundtrip tests parametrized with paths from another machine).pytest tests/test_vlasov1d/test_units_boundary.py— dimensional round-trip vs CODATA/NRLpytest tests/test_vlasov1d/test_fp_momentum_conservation.py— Dougherty conservation with n(x) ≠ 1vfp1dtest_kappa_ehstill passes (guards the untouchedvth_norm())