From 86341b9e32a0f6d7f37c2bcebd3f98af1c628246 Mon Sep 17 00:00:00 2001 From: Phil Travis Date: Fri, 17 Jul 2026 12:58:09 -0700 Subject: [PATCH 01/11] fix: adopt v0 = sqrt(T0/m_e) in electron_debye_normalization 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 --- adept/normalization.py | 14 ++-- tests/test_vlasov1d/test_units_boundary.py | 74 ++++++++++++++++++++++ 2 files changed, 83 insertions(+), 5 deletions(-) create mode 100644 tests/test_vlasov1d/test_units_boundary.py diff --git a/adept/normalization.py b/adept/normalization.py index aa9f7714..ae23bbc2 100644 --- a/adept/normalization.py +++ b/adept/normalization.py @@ -34,7 +34,7 @@ class PlasmaNormalization: def logLambda_ee(self) -> float: n0_cc = self.n0.to("1/cc").magnitude T0_eV = self.T0.to("eV").magnitude - logLambda_ee = 23.5 - jnp.log(n0_cc**0.5 / T0_eV**-1.25) + logLambda_ee = 23.5 - jnp.log(n0_cc**0.5 * T0_eV**-1.25) logLambda_ee -= (1e-5 + (jnp.log(T0_eV) - 2) ** 2.0 / 16) ** 0.5 return logLambda_ee @@ -78,9 +78,13 @@ def electron_debye_normalization(n0_str, T0_str): """ Returns the electron thermal normalization for the given density and temperature. Unit quantities are: - - Debye length - - Electron thermal velocity - - Langmuir oscillation frequency + - L0 = Debye length, sqrt(eps0 T0 / (n0 e^2)) + - v0 = electron thermal velocity sqrt(T0/m_e) (RMS / standard-deviation + convention: a T=1 Maxwellian is exp(-v^2/2) with unit variance) + - tau = 1/wp0 (inverse Langmuir oscillation frequency) + + Under this convention a code wavenumber is k*lambda_De and the Bohm-Gross + dispersion reads w^2 = 1 + 3 k^2. """ n0 = UREG.Quantity(n0_str) T0 = UREG.Quantity(T0_str) @@ -88,7 +92,7 @@ def electron_debye_normalization(n0_str, T0_str): wp0 = ((n0 * UREG.e**2.0 / (UREG.m_e * UREG.epsilon_0)) ** 0.5).to("rad/s") tau = 1 / wp0 - v0 = ((2.0 * T0 / UREG.m_e) ** 0.5).to("m/s") + v0 = ((T0 / UREG.m_e) ** 0.5).to("m/s") x0 = (v0 / wp0).to("nm") return PlasmaNormalization(m0=UREG.m_e, q0=UREG.e, n0=n0, T0=T0, L0=x0, v0=v0, tau=tau) diff --git a/tests/test_vlasov1d/test_units_boundary.py b/tests/test_vlasov1d/test_units_boundary.py new file mode 100644 index 00000000..e7af1c38 --- /dev/null +++ b/tests/test_vlasov1d/test_units_boundary.py @@ -0,0 +1,74 @@ +# Copyright (c) Ergodic LLC 2026 +# research@ergodic.io +""" +Units-boundary regression tests for the Vlasov-1D normalization layer. + +Every physics test in this suite works in dimensionless code units, which +validates the engine but not the dimensional dictionary: numeric inputs bypass +``normalize()`` entirely. These tests cross the units boundary — dimensional +input in, dimensional quantity out — and compare against *independent* physical +references (CODATA constants via scipy, and the NRL formulary), so a convention +error in ``normalization.py`` cannot hide. + +The engine convention under test: v0 = sqrt(T0/m_e) (RMS / standard-deviation +thermal speed), L0 = v0/wp0 = lambda_De, Maxwellian exp(-v^2/2) at T=1. +""" + +import numpy as np +from scipy import constants as csts + +from adept.normalization import electron_debye_normalization, normalize + +N0_STR, T0_STR = "1.5e21/cc", "2000eV" +N0_SI = 1.5e21 * 1e6 # 1/m^3 +T0_J = 2000.0 * csts.e + + +def test_debye_normalization_against_codata(): + """v0, L0, tau, and c_hat must match hand-computed CODATA values.""" + norm = electron_debye_normalization(N0_STR, T0_STR) + + v_th = np.sqrt(T0_J / csts.m_e) # sigma convention: sqrt(T/m) + wp0 = np.sqrt(N0_SI * csts.e**2 / (csts.epsilon_0 * csts.m_e)) + lambda_de = v_th / wp0 + + np.testing.assert_allclose(norm.v0.to("m/s").magnitude, v_th, rtol=1e-6) + np.testing.assert_allclose(norm.L0.to("m").magnitude, lambda_de, rtol=1e-6) + np.testing.assert_allclose(norm.tau.to("s").magnitude, 1.0 / wp0, rtol=1e-6) + np.testing.assert_allclose(norm.speed_of_light_norm(), csts.c / v_th, rtol=1e-6) + + +def test_debye_length_against_nrl_formulary(): + """lambda_De = 7.43e2 sqrt(T_eV/n_cc) cm (NRL formulary) — a truly external reference.""" + norm = electron_debye_normalization(N0_STR, T0_STR) + lambda_de_nrl_m = 7.43e2 * np.sqrt(2000.0 / 1.5e21) * 1e-2 + np.testing.assert_allclose(norm.L0.to("m").magnitude, lambda_de_nrl_m, rtol=1e-3) + + +def test_dimensional_string_round_trip(): + """String inputs must convert with L0 = lambda_De (not sqrt(2) lambda_De).""" + norm = electron_debye_normalization(N0_STR, T0_STR) + + v_th = np.sqrt(T0_J / csts.m_e) + wp0 = np.sqrt(N0_SI * csts.e**2 / (csts.epsilon_0 * csts.m_e)) + lambda_de = v_th / wp0 + + # x: 100 um -> 100e-6 / lambda_De code units + np.testing.assert_allclose(normalize("100um", norm, dim="x"), 100e-6 / lambda_de, rtol=1e-6) + # k: 1/um -> lambda_De / 1e-6 code units (k lambda_De) + np.testing.assert_allclose(normalize("1/um", norm, dim="k"), lambda_de / 1e-6, rtol=1e-6) + # t: 1 ps -> wp0 * 1e-12 + np.testing.assert_allclose(normalize("1ps", norm, dim="t"), wp0 * 1e-12, rtol=1e-6) + # numeric inputs pass through untouched + assert normalize(3.25, norm, dim="x") == 3.25 + + +def test_logged_collision_frequency_is_physical(): + """logLambda and nuee must be positive and match the NRL expression.""" + norm = electron_debye_normalization(N0_STR, T0_STR) + log_lambda = float(norm.logLambda_ee()) + # NRL: 23.5 - ln(n^1/2 T^-5/4) - [1e-5 + (ln T - 2)^2/16]^1/2 + expected = 23.5 - np.log(np.sqrt(1.5e21) * 2000.0**-1.25) - np.sqrt(1e-5 + (np.log(2000.0) - 2.0) ** 2 / 16.0) + np.testing.assert_allclose(log_lambda, expected, rtol=1e-10) + assert log_lambda > 0 + assert norm.approximate_ee_collision_frequency().to("Hz").magnitude > 0 From 2a9d687e9406a28b7ef3e904b9f43e5a29e4f0e5 Mon Sep 17 00:00:00 2001 From: Phil Travis Date: Fri, 17 Jul 2026 12:58:31 -0700 Subject: [PATCH 02/11] fix(tf1d): remove private sqrt(2) duplicate and logLambda sign error 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 --- adept/_tf1d/modules.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/adept/_tf1d/modules.py b/adept/_tf1d/modules.py index bdd53b4e..9816785d 100644 --- a/adept/_tf1d/modules.py +++ b/adept/_tf1d/modules.py @@ -59,7 +59,7 @@ def write_units(self): wp0 = np.sqrt(n0 * self.ureg.e**2.0 / (self.ureg.m_e * self.ureg.epsilon_0)).to("rad/s") tp0 = (1 / wp0).to("fs") - v0 = np.sqrt(2.0 * T0 / self.ureg.m_e).to("m/s") + v0 = np.sqrt(T0 / self.ureg.m_e).to("m/s") x0 = (v0 / wp0).to("nm") c_light = _Q(1.0 * self.ureg.c).to("m/s") / v0 beta = (v0 / self.ureg.c).to("dimensionless") @@ -72,7 +72,7 @@ def write_units(self): sim_duration = (self.cfg["grid"]["tmax"] * tp0).to("ps") # collisions - logLambda_ee = 23.5 - np.log(n0.magnitude**0.5 / T0.magnitude**-1.25) + logLambda_ee = 23.5 - np.log(n0.magnitude**0.5 * T0.magnitude**-1.25) logLambda_ee -= (1e-5 + (np.log(T0.magnitude) - 2) ** 2.0 / 16) ** 0.5 nuee = _Q(2.91e-6 * n0.magnitude * logLambda_ee / T0.magnitude**1.5, "Hz") nuee_norm = nuee / wp0 From c5563478a68953b059209141a6eb9de06ba9ff80 Mon Sep 17 00:00:00 2001 From: Phil Travis Date: Fri, 17 Jul 2026 12:58:32 -0700 Subject: [PATCH 03/11] fix(vlasov1d): Dougherty vbar missing 1/n; Krook target uses species 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 --- adept/_vlasov1d/modules.py | 6 ++ .../solvers/pushers/fokker_planck.py | 10 ++- .../test_fp_momentum_conservation.py | 83 +++++++++++++++++++ 3 files changed, 96 insertions(+), 3 deletions(-) create mode 100644 tests/test_vlasov1d/test_fp_momentum_conservation.py diff --git a/adept/_vlasov1d/modules.py b/adept/_vlasov1d/modules.py index fe2c52d3..40dcd090 100644 --- a/adept/_vlasov1d/modules.py +++ b/adept/_vlasov1d/modules.py @@ -117,6 +117,9 @@ def write_units(self) -> dict: sim_duration = (grid.tmax * norm.tau).to("ps") nu_ee = norm.approximate_ee_collision_frequency() + # e-e collision rate in code units (1/wp0). This is what the FP/Krook + # `baseline` rates should be compared against. + nu_ee_norm = (nu_ee * norm.tau).to("").magnitude beta = 1.0 / norm.speed_of_light_norm() @@ -130,6 +133,7 @@ def write_units(self) -> dict: "beta": beta, "x0": norm.L0.to("nm"), "nuee": nu_ee.to("Hz"), + "nuee_norm": nu_ee_norm, "logLambda_ee": norm.logLambda_ee(), "box_length": box_length, "box_width": box_width, @@ -235,10 +239,12 @@ def get_solver_quantities(self) -> dict: cfg_grid["species_grids"][species_name]["one_over_kvr"] = jnp.array(one_over_kvr) # Build species parameters (charge, mass, charge-to-mass ratio) + # T0 is the bulk (first-listed component) temperature, used e.g. by the Krook target cfg_grid["species_params"][species_name] = { "charge": species_cfg.charge, "mass": species_cfg.mass, "charge_to_mass": species_cfg.charge / species_cfg.mass, + "T0": self.simulation.species_distributions[species_name][0].T0, } cfg_grid["n_prof_total"] = n_prof_total diff --git a/adept/_vlasov1d/solvers/pushers/fokker_planck.py b/adept/_vlasov1d/solvers/pushers/fokker_planck.py index 23166bfd..dfb2a5b1 100644 --- a/adept/_vlasov1d/solvers/pushers/fokker_planck.py +++ b/adept/_vlasov1d/solvers/pushers/fokker_planck.py @@ -76,9 +76,9 @@ def compute_vbar(self, f: Array) -> Array: f: Distribution function, shape (..., nv) Returns: - vbar: Mean velocity , shape (...) + vbar: Mean velocity = ∫v f dv / ∫f dv, shape (...) """ - return jnp.sum(f * self.v, axis=-1) * self.dv + return jnp.sum(f * self.v, axis=-1) / jnp.sum(f, axis=-1) def compute_D(self, f: Array, beta: Array) -> Array: """ @@ -242,7 +242,11 @@ def __init__(self, cfg: Mapping[str, Any]): self.cfg = cfg v = cfg["grid"]["species_grids"]["electron"]["v"] dv = cfg["grid"]["species_grids"]["electron"]["dv"] - f_mx = np.exp(-(v[None, :] ** 2.0) / 2.0) + params = cfg["grid"].get("species_params", {}).get("electron", {}) + T0 = params.get("T0", 1.0) + mass = params.get("mass", 1.0) + # Maxwellian with variance T0/m (the species' bulk thermal width) + f_mx = np.exp(-(v[None, :] ** 2.0) / (2.0 * T0 / mass)) self.f_mx = f_mx / np.sum(f_mx, axis=1)[:, None] / dv self.dv = dv diff --git a/tests/test_vlasov1d/test_fp_momentum_conservation.py b/tests/test_vlasov1d/test_fp_momentum_conservation.py new file mode 100644 index 00000000..a50d9f66 --- /dev/null +++ b/tests/test_vlasov1d/test_fp_momentum_conservation.py @@ -0,0 +1,83 @@ +# Copyright (c) Ergodic LLC 2026 +# research@ergodic.io +""" +Momentum conservation of the Dougherty operator where n(x) != 1. + +The Dougherty drag must center on the *mean velocity* u(x) = ∫v f dv / n(x). +A drag centered on the raw first moment ∫v f dv = n*u (the historical bug) +breaks momentum conservation exactly where the density deviates from 1, e.g. +for bump-on-tail or large-amplitude density perturbations with collisions on. + +This test applies the collision operator in isolation (no advection, no +fields) to a drifting Maxwellian with a strong density modulation and checks +per-cell conservation of density, momentum, and energy. +""" + +import jax + +jax.config.update("jax_enable_x64", True) + +import numpy as np +import pytest +from jax import numpy as jnp + +from adept._vlasov1d.solvers.pushers.fokker_planck import Collisions + + +def _make_cfg(nx, nv, vmax, fp_type): + dv = 2.0 * vmax / nv + v = jnp.linspace(-vmax + dv / 2.0, vmax - dv / 2.0, nv) + return { + "grid": { + "species_grids": {"electron": {"v": v, "dv": dv, "nv": nv, "vmax": vmax}}, + "species_params": {"electron": {"charge": -1.0, "mass": 1.0, "charge_to_mass": -1.0, "T0": 1.0}}, + }, + "terms": { + "fokker_planck": {"is_on": True, "type": fp_type}, + "krook": {"is_on": False}, + }, + } + + +# Energy tolerance is per-scheme: Chang-Cooper preserves the discrete equilibrium +# to ~1e-6, plain central differencing only to ~1e-3 at this resolution. +@pytest.mark.parametrize("fp_type,energy_rtol", [("dougherty", 5e-3), ("chang_cooper_dougherty", 1e-5)]) +def test_dougherty_momentum_conservation_with_density_perturbation(fp_type, energy_rtol): + nx, nv, vmax = 16, 512, 6.4 + cfg = _make_cfg(nx, nv, vmax, fp_type) + v = np.array(cfg["grid"]["species_grids"]["electron"]["v"]) + dv = cfg["grid"]["species_grids"]["electron"]["dv"] + + # Strong density modulation and a finite drift: exactly the regime where + # centering the drag on n*u instead of u destroys momentum conservation. + x = np.linspace(0, 2 * np.pi, nx, endpoint=False) + n_prof = 1.0 + 0.5 * np.sin(x) + u0 = 0.5 + f = n_prof[:, None] * np.exp(-((v[None, :] - u0) ** 2) / 2.0) + f = f / (np.sum(np.exp(-((v - u0) ** 2) / 2.0)) * dv) + + f = jnp.array(f) + collisions = Collisions(cfg) + + nu = jnp.ones(nx) + dt = 0.1 + f_out = f + for _ in range(50): + f_out = collisions(nu, None, f_out, dt) + + n0_ = np.sum(np.array(f), axis=1) * dv + n1 = np.sum(np.array(f_out), axis=1) * dv + p0 = np.sum(np.array(f) * v[None, :], axis=1) * dv + p1 = np.sum(np.array(f_out) * v[None, :], axis=1) * dv + e0 = np.sum(np.array(f) * v[None, :] ** 2, axis=1) * dv + e1 = np.sum(np.array(f_out) * v[None, :] ** 2, axis=1) * dv + + # The historical n*u-centered drag gives O(50%) momentum errors here; + # the corrected operator conserves to discretization level (~1e-7). + np.testing.assert_allclose(n1, n0_, rtol=1e-10, err_msg="density not conserved") + np.testing.assert_allclose(p1, p0, rtol=1e-6, err_msg="momentum not conserved where n(x) != 1") + np.testing.assert_allclose(e1, e0, rtol=energy_rtol, err_msg="energy not conserved") + + # And the distribution must still be centered on u0, not n(x)*u0 + u_final = p1 / n1 + np.testing.assert_allclose(u_final, u0, atol=1e-6) From 1ea7244af6d918c725ae0df6ffac5dd2c84b15d4 Mon Sep 17 00:00:00 2001 From: Phil Travis Date: Fri, 17 Jul 2026 12:58:32 -0700 Subject: [PATCH 04/11] fix(vlasov1d): center storage moments on u not n*u; entropy sign; energy 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 --- adept/_vlasov1d/storage.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/adept/_vlasov1d/storage.py b/adept/_vlasov1d/storage.py index 1075ed48..6899e049 100644 --- a/adept/_vlasov1d/storage.py +++ b/adept/_vlasov1d/storage.py @@ -137,11 +137,12 @@ def _calc_moment_(inp, _dv=dv): f = y[species_name] species_moments = {} species_moments["n"] = _calc_moment_(f) - species_moments["v"] = _calc_moment_(f * v[None, :]) + species_moments["j"] = _calc_moment_(f * v[None, :]) + species_moments["v"] = species_moments["j"] / species_moments["n"] v_m_vbar = v[None, :] - species_moments["v"][:, None] species_moments["p"] = _calc_moment_(f * v_m_vbar**2.0) species_moments["q"] = _calc_moment_(f * v_m_vbar**3.0) - species_moments["-flogf"] = _calc_moment_(f * jnp.log(jnp.abs(f))) + species_moments["-flogf"] = _calc_moment_(-jnp.abs(f) * jnp.log(jnp.abs(f))) species_moments["f^2"] = _calc_moment_(f * f) result[species_name] = species_moments @@ -292,9 +293,11 @@ def save(t, y, args): scalars = {} # Compute scalars for each species + mean_kinetic_energy = 0.0 for species_name in species_names: v = species_grids[species_name]["v"][None, :] dv = species_grids[species_name]["dv"] + mass = cfg["grid"]["species_params"][species_name]["mass"] def _calc_mean_moment_(inp, _dv=dv): return jnp.mean(jnp.sum(inp, axis=1) * _dv) @@ -306,12 +309,19 @@ def _calc_mean_moment_(inp, _dv=dv): scalars[f"mean_q_{species_name}"] = _calc_mean_moment_(f * v**3.0) scalars[f"mean_-flogf_{species_name}"] = _calc_mean_moment_(-jnp.log(jnp.abs(f)) * jnp.abs(f)) scalars[f"mean_f2_{species_name}"] = _calc_mean_moment_(f * f) + mean_kinetic_energy += 0.5 * mass * scalars[f"mean_P_{species_name}"] # Shared field scalars (not species-specific) scalars["mean_de2"] = jnp.mean(y["de"] ** 2.0) scalars["mean_e2"] = jnp.mean(y["e"] ** 2.0) scalars["mean_pond"] = jnp.mean(-0.5 * jnp.gradient(y["a"] ** 2.0, cfg["grid"]["dx"])[1:-1]) + # Energy conservation monitor (electrostatic): x-averaged kinetic + E-field + # energy density. Transverse (EM) field energy is not included. + scalars["mean_kinetic_energy"] = mean_kinetic_energy + scalars["mean_field_energy"] = 0.5 * scalars["mean_e2"] + scalars["mean_total_energy"] = mean_kinetic_energy + 0.5 * scalars["mean_e2"] + return scalars return save From 070fb044cf994da3b1b77d562a597e5821a7bf5f Mon Sep 17 00:00:00 2001 From: Phil Travis Date: Fri, 17 Jul 2026 12:58:32 -0700 Subject: [PATCH 05/11] fix(vlasov1d): dx and semi-Lagrangian period ignored xmin 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 --- adept/_vlasov1d/grid.py | 2 +- adept/_vlasov1d/solvers/pushers/vlasov.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/adept/_vlasov1d/grid.py b/adept/_vlasov1d/grid.py index 91b6d09b..b16db67e 100644 --- a/adept/_vlasov1d/grid.py +++ b/adept/_vlasov1d/grid.py @@ -52,7 +52,7 @@ def __init__( self.tmin = tmin # Compute dx - self.dx = xmax / nx + self.dx = (xmax - xmin) / nx # Override dt for EM wave stability if needed if should_override_dt_for_em_waves: diff --git a/adept/_vlasov1d/solvers/pushers/vlasov.py b/adept/_vlasov1d/solvers/pushers/vlasov.py index 200da317..4a51ec4c 100644 --- a/adept/_vlasov1d/solvers/pushers/vlasov.py +++ b/adept/_vlasov1d/solvers/pushers/vlasov.py @@ -28,7 +28,7 @@ def __init__(self, cfg, interp_e): """Create interpolation helpers for x-advection and v-advection steps.""" self.x = cfg["grid"]["x"] self.v = cfg["grid"]["v"] - self.f_interp = partial(interp2d, x=self.x, y=self.v, period=cfg["grid"]["xmax"]) + self.f_interp = partial(interp2d, x=self.x, y=self.v, period=cfg["grid"]["xmax"] - cfg["grid"]["xmin"]) self.dt = cfg["grid"]["dt"] self.dummy_x = jnp.ones_like(self.x) self.dummy_v = jnp.ones_like(self.v)[None, :] From 0b626bb2529e94a74459425be654f69d20f320ca Mon Sep 17 00:00:00 2001 From: Phil Travis Date: Fri, 17 Jul 2026 12:58:32 -0700 Subject: [PATCH 06/11] chore(vlasov1d): remove dead c_light knob; fix AmpereSolver docstring 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 --- adept/_vlasov1d/datamodel.py | 1 - adept/_vlasov1d/solvers/pushers/field.py | 2 +- configs/vlasov-1d/wavepacket.yaml | 1 - 3 files changed, 1 insertion(+), 3 deletions(-) diff --git a/adept/_vlasov1d/datamodel.py b/adept/_vlasov1d/datamodel.py index 2e96c1ca..811df767 100644 --- a/adept/_vlasov1d/datamodel.py +++ b/adept/_vlasov1d/datamodel.py @@ -70,7 +70,6 @@ class GridConfig(BaseModel): nv: int | None = None vmax: float | None = None parallel: tuple[str, ...] | bool = False - c_light: float | None = None @field_validator("parallel", mode="before") @classmethod diff --git a/adept/_vlasov1d/solvers/pushers/field.py b/adept/_vlasov1d/solvers/pushers/field.py index 76a999af..ff73262c 100644 --- a/adept/_vlasov1d/solvers/pushers/field.py +++ b/adept/_vlasov1d/solvers/pushers/field.py @@ -227,7 +227,7 @@ def __call__(self, f_dict: dict, prev_ex: jnp.ndarray, dt: jnp.float64): class AmpereSolver: """Ampere solver using current density to evolve electric field. - Solves: a single Euler step of ∂E/∂t = -j, where j = Σ_s (q_s/m_s) ∫v f_s dv + Solves: a single Euler step of ∂E/∂t = -j, where j = Σ_s q_s ∫v f_s dv is the total current density from all species. """ diff --git a/configs/vlasov-1d/wavepacket.yaml b/configs/vlasov-1d/wavepacket.yaml index 3770512f..5ceba970 100644 --- a/configs/vlasov-1d/wavepacket.yaml +++ b/configs/vlasov-1d/wavepacket.yaml @@ -35,7 +35,6 @@ drivers: width: 900.0 ey: {} grid: - c_light: 10 dt: 0.1 nv: 6144 nx: 4096 From 438d4b4050c05041a3988d8e375a452ec759c2fa Mon Sep 17 00:00:00 2001 From: Phil Travis Date: Fri, 17 Jul 2026 12:58:54 -0700 Subject: [PATCH 07/11] fix(vlasov1d): retune srs-debug-small numeric driver to corrected c_hat 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 --- configs/vlasov-1d/srs-debug-small.yaml | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/configs/vlasov-1d/srs-debug-small.yaml b/configs/vlasov-1d/srs-debug-small.yaml index c0bfa5ec..5ac73e31 100644 --- a/configs/vlasov-1d/srs-debug-small.yaml +++ b/configs/vlasov-1d/srs-debug-small.yaml @@ -48,13 +48,18 @@ mlflow: experiment: vlasov1d-srs run: debug-small +# The ey driver pair is a backscatter-SRS matched triad in code units +# (v0 = sqrt(T0/m_e), k in 1/lambda_De, w in wp0, c_hat = c/v0 = 15.984 at 2000 eV): +# pump: w0 = 2.79, k0 = +sqrt(w0^2 - 1)/c_hat = +0.162949 (rightgoing) +# seed: w0 = 1.700942, k0 = -0.086080 (leftgoing) +# EPW: k lambda_De = 0.249029, w = 1.089058 = w_pump - w_seed drivers: ex: {} ey: '0': params: a0: 1.0 - k0: 1.0 + k0: 0.162949 w0: 2.79 dw0: 0. envelope: @@ -69,8 +74,8 @@ drivers: '1': params: a0: 1.e-6 - k0: 1.0 - w0: 1.63 + k0: -0.086080 + w0: 1.700942 dw0: 0. envelope: time: From cc7fd581141931891ac8b7f2b7ec042b4d83280c Mon Sep 17 00:00:00 2001 From: Phil Travis Date: Fri, 17 Jul 2026 12:58:54 -0700 Subject: [PATCH 08/11] docs(vlasov1d): state the normalization convention; document unit traps - 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 / = 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 --- adept/_vlasov1d/helpers.py | 5 ++++- docs/source/solvers/vlasov1d/config.md | 30 ++++++++++++++++++++++---- 2 files changed, 30 insertions(+), 5 deletions(-) diff --git a/adept/_vlasov1d/helpers.py b/adept/_vlasov1d/helpers.py index 553c8311..5a37afcf 100644 --- a/adept/_vlasov1d/helpers.py +++ b/adept/_vlasov1d/helpers.py @@ -68,7 +68,10 @@ def _initialize_supergaussian_distribution_( # Thermal velocity: v_t = sqrt(T/m) v_thermal = np.sqrt(T0 / mass) - # Alpha factor for supergaussian normalization + # Alpha factor for supergaussian normalization. This fixes the moment ratio + # / = 3*T0/mass for every order m; the realized *variance* equals + # T0/mass only at m=2 (Maxwellian). For m>2 flat-tops the second-moment + # temperature exceeds T0 (x1.24 at m=3, x1.37 at m=4). See docs config.md. alpha = np.sqrt(3.0 * gamma_3_over_m(supergaussian_order) / gamma_5_over_m(supergaussian_order)) single_dist = -(np.power(np.abs((vax[None, :] - v0) / (alpha * v_thermal)), supergaussian_order)) diff --git a/docs/source/solvers/vlasov1d/config.md b/docs/source/solvers/vlasov1d/config.md index 66e878f7..f092c65e 100644 --- a/docs/source/solvers/vlasov1d/config.md +++ b/docs/source/solvers/vlasov1d/config.md @@ -48,6 +48,28 @@ units: normalizing_density: 1.5e21/cc ``` +### Normalization convention + +All Vlasov-1D quantities are normalized using the following unit system, built from +`normalizing_density` ($n_0$) and `normalizing_temperature` ($T_0$): + +| Unit | Definition | Meaning | +|------|-----------|---------| +| time | $\tau = 1/\omega_{p0}$, $\omega_{p0} = \sqrt{n_0 e^2/(\epsilon_0 m_e)}$ | inverse plasma frequency | +| velocity | $v_0 = \sqrt{T_0/m_e}$ | electron thermal speed (RMS / standard-deviation convention) | +| length | $L_0 = v_0/\omega_{p0} = \lambda_{De}$ | electron Debye length | +| wavenumber | $1/L_0$ | a code wavenumber is $k \lambda_{De}$ | + +Consequences of the $v_0 = \sqrt{T_0/m_e}$ (σ) convention: + +- A Maxwellian at code temperature $T = 1$ is $f \propto e^{-v^2/2}$ (unit variance). +- The Bohm–Gross dispersion in code units is $\omega^2 = 1 + 3 k^2$. +- The normalized speed of light is $\hat c = c/\sqrt{T_0/m_e}$ (e.g. $\hat c = 15.98$ at 2000 eV). + +Dimensional inputs (strings with units, e.g. `xmax: 100um`) are converted with these +units; plain numeric inputs are taken to already be in code units and pass through +unchanged. + ## density Species and density configuration. You can define multiple "species" by using keys prefixed with `species-`. @@ -65,9 +87,9 @@ Each species is defined with a key starting with `species-` (e.g., `species-back | `noise_seed` | int | Random seed for noise initialization | | `noise_type` | string | `"gaussian"` or `"uniform"` | | `noise_val` | float | Amplitude of noise | -| `v0` | float | Drift velocity (normalized) | -| `T0` | float | Temperature (normalized to `normalizing_temperature`) | -| `m` | float | Exponent for super-Gaussian distribution. `2.0` is Maxwellian | +| `v0` | float | Drift velocity in code units of $\sqrt{T_0/m_e}$ (thermal-σ units). Numeric only — dimensional strings are not supported here | +| `T0` | float | Temperature in units of `normalizing_temperature`. The initialized distribution has velocity variance `T0/mass`. Numeric only | +| `m` | float | Exponent for super-Gaussian distribution $f \propto \exp[-\|v/(\alpha v_{th})\|^m]$. `2.0` is Maxwellian. Note: $\alpha$ is chosen to fix the moment ratio $\langle v^4\rangle/\langle v^2\rangle = 3\,T_0/\mathrm{mass}$ for all $m$; the *variance* equals `T0/mass` only at `m: 2`. For flat-top distributions (`m > 2`) the second-moment temperature diagnostic will read higher than `T0` (e.g. ×1.24 at `m: 3`, ×1.37 at `m: 4`) | | `basis` | string | Spatial profile type (see below) | #### Basis Types @@ -531,7 +553,7 @@ Both `fokker_planck` and `krook` use profile objects to define spatiotemporal va | Field | Type | Description | |-------|------|-------------| -| `baseline` | float | Base collision frequency | +| `baseline` | float | Base collision frequency in code units of $\omega_{p0}$ (i.e. $\hat\nu = \nu[\mathrm{s}^{-1}]/\omega_{p0}$, a true rate — no $2\pi$). For reference, the NRL e–e rate in these units is logged as `nuee_norm` in `units.yaml`. Note the Dougherty/Lenard–Bernstein $\nu$ agrees with the NRL $\nu_{ee}$ only up to an O(1) factor | | `bump_or_trough` | string | `"bump"` or `"trough"` | | `center` | float | Profile center | | `rise` | float | Transition steepness | From f9dfb4c085d02d35e64bd98f7fc85f75cc04d4af Mon Sep 17 00:00:00 2001 From: Phil Travis Date: Fri, 17 Jul 2026 12:58:55 -0700 Subject: [PATCH 09/11] test(vlasov1d): regenerate config regression fixtures 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 --- ...okker_planck_conservation_array_config.yml | 18 ++++++++++-------- ...ker_planck_conservation_derived_config.yml | 17 +++++++++-------- .../fokker_planck_conservation_units.yml | 15 ++++++++------- ...multispecies_ion_acoustic_array_config.yml | 19 +++++++++++-------- ...ltispecies_ion_acoustic_derived_config.yml | 17 +++++++++-------- .../multispecies_ion_acoustic_units.yml | 15 ++++++++------- .../resonance_array_config.yml | 18 ++++++++++-------- .../resonance_derived_config.yml | 17 +++++++++-------- .../resonance_units.yml | 15 ++++++++------- 9 files changed, 82 insertions(+), 69 deletions(-) diff --git a/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_array_config.yml b/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_array_config.yml index 4134829c..9539860e 100644 --- a/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_array_config.yml +++ b/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_array_config.yml @@ -38,7 +38,7 @@ drivers: w0: 1.0 ey: {} grid: - beta: 0.088474881879774 + beta: 0.062561188941867 dt: 0.1 dx: 0.785 ion_charge: @@ -16507,6 +16507,7 @@ grid: vmax: 6.4 species_params: electron: + T0: 1.0 charge: -1.0 charge_to_mass: -1.0 mass: 1.0 @@ -17723,17 +17724,18 @@ terms: units: derived: T0: 2000 electron_volt - beta: 0.088474881879774 - box_length: 0.07623654008768617 micron + beta: 0.062561188941867 + box_length: 0.05390737447020297 micron box_width: inf - c_light: 11.302642950785 - logLambda_ee: -11.781233290653 + c_light: 15.984350951661 + logLambda_ee: 7.221022858202 n0: 1.5e+21 / cubic_centimeter - nuee: -574949910190.1326 hertz + nuee: 352401683370.27344 hertz + nuee_norm: 0.00016128753860756 sim_duration: 0.004622577634581562 picosecond tp0: 0.45768095391896646 femtosecond - v0: 26524102.30999721 meter / second + v0: 18755372.608284798 meter / second wp0: 2184928150138955.5 radian / second - x0: 12.139576447083783 nanometer + x0: 8.58397682646544 nanometer normalizing_density: 1.5e21/cc normalizing_temperature: 2000eV diff --git a/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_derived_config.yml b/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_derived_config.yml index ebca6209..f83b3fe6 100644 --- a/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_derived_config.yml +++ b/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_derived_config.yml @@ -38,7 +38,7 @@ drivers: w0: 1.0 ey: {} grid: - beta: 0.088474881879774 + beta: 0.062561188941867 dt: 0.1 dx: 0.785 max_steps: 105 @@ -111,17 +111,18 @@ terms: units: derived: T0: 2000 electron_volt - beta: 0.088474881879774 - box_length: 0.07623654008768617 micron + beta: 0.062561188941867 + box_length: 0.05390737447020297 micron box_width: inf - c_light: 11.302642950785 - logLambda_ee: -11.781233290653 + c_light: 15.984350951661 + logLambda_ee: 7.221022858202 n0: 1.5e+21 / cubic_centimeter - nuee: -574949910190.1326 hertz + nuee: 352401683370.27344 hertz + nuee_norm: 0.00016128753860756 sim_duration: 0.004622577634581562 picosecond tp0: 0.45768095391896646 femtosecond - v0: 26524102.30999721 meter / second + v0: 18755372.608284798 meter / second wp0: 2184928150138955.5 radian / second - x0: 12.139576447083783 nanometer + x0: 8.58397682646544 nanometer normalizing_density: 1.5e21/cc normalizing_temperature: 2000eV diff --git a/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_units.yml b/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_units.yml index dca5edde..f1e7ae4b 100644 --- a/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_units.yml +++ b/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_units.yml @@ -1,13 +1,14 @@ T0: 2000 electron_volt -beta: 0.088474881879774 -box_length: 0.07623654008768617 micron +beta: 0.062561188941867 +box_length: 0.05390737447020297 micron box_width: inf -c_light: 11.302642950785 -logLambda_ee: -11.781233290653 +c_light: 15.984350951661 +logLambda_ee: 7.221022858202 n0: 1.5e+21 / cubic_centimeter -nuee: -574949910190.1326 hertz +nuee: 352401683370.27344 hertz +nuee_norm: 0.00016128753860756 sim_duration: 0.004622577634581562 picosecond tp0: 0.45768095391896646 femtosecond -v0: 26524102.30999721 meter / second +v0: 18755372.608284798 meter / second wp0: 2184928150138955.5 radian / second -x0: 12.139576447083783 nanometer +x0: 8.58397682646544 nanometer diff --git a/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_array_config.yml b/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_array_config.yml index 927d9760..e899d7fa 100644 --- a/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_array_config.yml +++ b/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_array_config.yml @@ -29,7 +29,7 @@ drivers: ex: {} ey: {} grid: - beta: 0.088474881879774 + beta: 0.062561188941867 dt: 0.5 dx: 9.81746875 ion_charge: @@ -28711,10 +28711,12 @@ grid: vmax: 0.005 species_params: electron: + T0: 1.0 charge: -1.0 charge_to_mass: -1.0 mass: 1.0 ion: + T0: 0.01 charge: 10.0 charge_to_mass: 0.00054466230936819 mass: 18360.0 @@ -38873,17 +38875,18 @@ terms: units: derived: T0: 2000 electron_volt - beta: 0.088474881879774 - box_length: 3.8137571970393944 micron + beta: 0.062561188941867 + box_length: 2.696733575825556 micron box_width: inf - c_light: 11.302642950785 - logLambda_ee: -11.781233290653 + c_light: 15.984350951661 + logLambda_ee: 7.221022858202 n0: 1.5e+21 / cubic_centimeter - nuee: -574949910190.1326 hertz + nuee: 352401683370.27344 hertz + nuee_norm: 0.00016128753860756 sim_duration: 2.288633610071792 picosecond tp0: 0.45768095391896646 femtosecond - v0: 26524102.30999721 meter / second + v0: 18755372.608284798 meter / second wp0: 2184928150138955.5 radian / second - x0: 12.139576447083783 nanometer + x0: 8.58397682646544 nanometer normalizing_density: 1.5e21/cc normalizing_temperature: 2000eV diff --git a/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_derived_config.yml b/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_derived_config.yml index e7c61532..9fe52c26 100644 --- a/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_derived_config.yml +++ b/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_derived_config.yml @@ -29,7 +29,7 @@ drivers: ex: {} ey: {} grid: - beta: 0.088474881879774 + beta: 0.062561188941867 dt: 0.5 dx: 9.81746875 max_steps: 10005 @@ -120,17 +120,18 @@ terms: units: derived: T0: 2000 electron_volt - beta: 0.088474881879774 - box_length: 3.8137571970393944 micron + beta: 0.062561188941867 + box_length: 2.696733575825556 micron box_width: inf - c_light: 11.302642950785 - logLambda_ee: -11.781233290653 + c_light: 15.984350951661 + logLambda_ee: 7.221022858202 n0: 1.5e+21 / cubic_centimeter - nuee: -574949910190.1326 hertz + nuee: 352401683370.27344 hertz + nuee_norm: 0.00016128753860756 sim_duration: 2.288633610071792 picosecond tp0: 0.45768095391896646 femtosecond - v0: 26524102.30999721 meter / second + v0: 18755372.608284798 meter / second wp0: 2184928150138955.5 radian / second - x0: 12.139576447083783 nanometer + x0: 8.58397682646544 nanometer normalizing_density: 1.5e21/cc normalizing_temperature: 2000eV diff --git a/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_units.yml b/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_units.yml index b05d7d96..e018602a 100644 --- a/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_units.yml +++ b/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_units.yml @@ -1,13 +1,14 @@ T0: 2000 electron_volt -beta: 0.088474881879774 -box_length: 3.8137571970393944 micron +beta: 0.062561188941867 +box_length: 2.696733575825556 micron box_width: inf -c_light: 11.302642950785 -logLambda_ee: -11.781233290653 +c_light: 15.984350951661 +logLambda_ee: 7.221022858202 n0: 1.5e+21 / cubic_centimeter -nuee: -574949910190.1326 hertz +nuee: 352401683370.27344 hertz +nuee_norm: 0.00016128753860756 sim_duration: 2.288633610071792 picosecond tp0: 0.45768095391896646 femtosecond -v0: 26524102.30999721 meter / second +v0: 18755372.608284798 meter / second wp0: 2184928150138955.5 radian / second -x0: 12.139576447083783 nanometer +x0: 8.58397682646544 nanometer diff --git a/tests/test_vlasov1d/test_config_regression/resonance_array_config.yml b/tests/test_vlasov1d/test_config_regression/resonance_array_config.yml index 8912206a..0924eec1 100644 --- a/tests/test_vlasov1d/test_config_regression/resonance_array_config.yml +++ b/tests/test_vlasov1d/test_config_regression/resonance_array_config.yml @@ -38,7 +38,7 @@ drivers: w0: 1.1598 ey: {} grid: - beta: 0.088474881879774 + beta: 0.062561188941867 dt: 0.25 dx: 0.654375 ion_charge: @@ -20747,6 +20747,7 @@ grid: vmax: 6.4 species_params: electron: + T0: 1.0 charge: -1.0 charge_to_mass: -1.0 mass: 1.0 @@ -23327,17 +23328,18 @@ terms: units: derived: T0: 2000 electron_volt - beta: 0.088474881879774 - box_length: 0.25420273080193445 micron + beta: 0.062561188941867 + box_length: 0.17974847474618633 micron box_width: inf - c_light: 11.302642950785 - logLambda_ee: -11.781233290653 + c_light: 15.984350951661 + logLambda_ee: 7.221022858202 n0: 1.5e+21 / cubic_centimeter - nuee: -574949910190.1326 hertz + nuee: 352401683370.27344 hertz + nuee_norm: 0.00016128753860756 sim_duration: 0.21980127811958367 picosecond tp0: 0.45768095391896646 femtosecond - v0: 26524102.30999721 meter / second + v0: 18755372.608284798 meter / second wp0: 2184928150138955.5 radian / second - x0: 12.139576447083783 nanometer + x0: 8.58397682646544 nanometer normalizing_density: 1.5e21/cc normalizing_temperature: 2000eV diff --git a/tests/test_vlasov1d/test_config_regression/resonance_derived_config.yml b/tests/test_vlasov1d/test_config_regression/resonance_derived_config.yml index 66794372..6e5e2484 100644 --- a/tests/test_vlasov1d/test_config_regression/resonance_derived_config.yml +++ b/tests/test_vlasov1d/test_config_regression/resonance_derived_config.yml @@ -38,7 +38,7 @@ drivers: w0: 1.1598 ey: {} grid: - beta: 0.088474881879774 + beta: 0.062561188941867 dt: 0.25 dx: 0.654375 max_steps: 1925 @@ -119,17 +119,18 @@ terms: units: derived: T0: 2000 electron_volt - beta: 0.088474881879774 - box_length: 0.25420273080193445 micron + beta: 0.062561188941867 + box_length: 0.17974847474618633 micron box_width: inf - c_light: 11.302642950785 - logLambda_ee: -11.781233290653 + c_light: 15.984350951661 + logLambda_ee: 7.221022858202 n0: 1.5e+21 / cubic_centimeter - nuee: -574949910190.1326 hertz + nuee: 352401683370.27344 hertz + nuee_norm: 0.00016128753860756 sim_duration: 0.21980127811958367 picosecond tp0: 0.45768095391896646 femtosecond - v0: 26524102.30999721 meter / second + v0: 18755372.608284798 meter / second wp0: 2184928150138955.5 radian / second - x0: 12.139576447083783 nanometer + x0: 8.58397682646544 nanometer normalizing_density: 1.5e21/cc normalizing_temperature: 2000eV diff --git a/tests/test_vlasov1d/test_config_regression/resonance_units.yml b/tests/test_vlasov1d/test_config_regression/resonance_units.yml index cc18e176..22e9776c 100644 --- a/tests/test_vlasov1d/test_config_regression/resonance_units.yml +++ b/tests/test_vlasov1d/test_config_regression/resonance_units.yml @@ -1,13 +1,14 @@ T0: 2000 electron_volt -beta: 0.088474881879774 -box_length: 0.25420273080193445 micron +beta: 0.062561188941867 +box_length: 0.17974847474618633 micron box_width: inf -c_light: 11.302642950785 -logLambda_ee: -11.781233290653 +c_light: 15.984350951661 +logLambda_ee: 7.221022858202 n0: 1.5e+21 / cubic_centimeter -nuee: -574949910190.1326 hertz +nuee: 352401683370.27344 hertz +nuee_norm: 0.00016128753860756 sim_duration: 0.21980127811958367 picosecond tp0: 0.45768095391896646 femtosecond -v0: 26524102.30999721 meter / second +v0: 18755372.608284798 meter / second wp0: 2184928150138955.5 radian / second -x0: 12.139576447083783 nanometer +x0: 8.58397682646544 nanometer From f5785ad28c9b08521b41a594e3221f31dcb4d84c Mon Sep 17 00:00:00 2001 From: Phil Travis Date: Mon, 20 Jul 2026 11:44:10 -0700 Subject: [PATCH 10/11] Added audit files --- .../2026-07-16_ADEPT_CONVENTIONS_AUDIT.md | 225 ++++++++++++++++ .../2026-07-16_VLASOV1D_CONVENTIONS_AUDIT.md | 249 ++++++++++++++++++ 2 files changed, 474 insertions(+) create mode 100644 docs/audits/2026-07-16_ADEPT_CONVENTIONS_AUDIT.md create mode 100644 docs/audits/2026-07-16_VLASOV1D_CONVENTIONS_AUDIT.md diff --git a/docs/audits/2026-07-16_ADEPT_CONVENTIONS_AUDIT.md b/docs/audits/2026-07-16_ADEPT_CONVENTIONS_AUDIT.md new file mode 100644 index 00000000..4b1ac623 --- /dev/null +++ b/docs/audits/2026-07-16_ADEPT_CONVENTIONS_AUDIT.md @@ -0,0 +1,225 @@ +# ADEPT Codebase-Wide Physics & Conventions Audit + +**Date:** 2026-07-16 +**Scope:** every solver module and shared physics file in `adept/` — `_vlasov2d`, `_pic1d`, `_tf1d`, `_lpse2d`, `_spectrax1d`, `_hermite_poisson_1d`, `vfp1d`, `osiris`, `vlasov1d2v`, plus the shared `normalization.py`, `electrostatic.py`, `functions.py`, `driftdiffusion.py`, `utils.py`, and all tests/configs. +**Companion document:** `VLASOV1D_CONVENTIONS_AUDIT.md` (the `_vlasov1d` deep audit that established the root cause). This document extends that audit to the rest of the codebase and answers: *which solvers are immune?* +**Method:** six parallel independent module audits, each cross-checked against source; load-bearing claims re-verified by hand. No source files modified. + +--- + +## 1. Executive summary + +**The √2 bug is confined to the `electron_debye_normalization` consumers — and even there, only to the unit-conversion/reporting layer.** The dimensionless *dynamics* of every solver in the codebase are internally self-consistent. The full damage inventory of `normalization.py:91` ($v_0 = \sqrt{2T_0/m_e}$, $L_0 = \sqrt2\lambda_{De}$): + +- `_vlasov1d`, `_vlasov2d`, `_pic1d`: logged physical units √2-wrong, physical temperature 2× the label, dimensional-string inputs √2-off, EM-branch $\hat c$ √2-small (details per module below). +- `_tf1d`: **independently re-implements the identical √2** in its own `write_units` (`_tf1d/modules.py:62`) — same symptom, separate code, needs its own fix. +- Everyone else: **immune** (see the table). + +**One critical amendment to the fix prescription in the vlasov1d audit:** `vth_norm()` (`normalization.py:48–50`) must **NOT** be redefined. Grep-verified, it has exactly four callers, all in `vfp1d/`, which is built self-consistently on the $v_{th} \equiv \sqrt{2T/m}$ (most-probable-speed) convention — its initializer carries a compensating factor (`vfp1d/helpers.py:111`: $\alpha = \sqrt{3\Gamma(3/m)/(2\Gamma(5/m))} = 1$ at $m{=}2$, vs. `_vlasov1d`'s $\sqrt2$) so the realized variance is exactly $T/m$. Dropping the 2 in `vth_norm()` would do nothing for `_vlasov1d` (which never calls it) and would **silently halve VFP-1D's initialized temperature** while leaving its transport coefficients at $T_0$. The safe fix touches only `electron_debye_normalization` lines 91–92. + +### 1.1 Immunity table (the answer to "which solvers are immune") + +"Bug" = the current √2 in `normalization.py:91`. "Fix" = changing `electron_debye_normalization` to $v_0=\sqrt{T_0/m_e}$ (only; `vth_norm()` untouched). + +| Solver | Working vth convention | Affected by bug NOW? | Affected by fix? | Notes | +|---|---|---|---|---| +| `_vlasov1d` | $\sqrt{T/m}$ (σ) | **YES** — logged units, physical T (2×), dimensional strings, EM $\hat c$; fixtures lock wrong values | Fix-safe for tests; must regenerate 3 fixtures | See companion doc | +| `_vlasov2d` | $\sqrt{T/m}$ (σ) | **YES (partial)** — same chain: logged units, physical T, EM $\hat c$; ES dynamics immune | **Fix-safe** — all 3 tests self-consistent, **no fixtures** | Cleanest of the affected: already fixed several vlasov1d bugs | +| `_pic1d` | $\sqrt{T/m}$ (σ) | **Logged units only** — dynamics immune (all inputs numeric) | **Fix-safe** — no fixtures, no test depends on it | Loader ≡ vlasov1d initializer convention | +| `_tf1d` | $\sqrt{T/m}$ (σ) | Not via `normalization.py` — but its **own private copy** of the √2 (`modules.py:62`) corrupts its logged units identically | Fix does **nothing** here; `modules.py:62` needs its own `2.0*T0 → T0` | Dynamics/tests fully decoupled from both | +| `_lpse2d` | $\sqrt{T/m}$ (σ), physical units | **IMMUNE** — never touches `normalization.py` | **IMMUNE** | Own ps/µm unit system; Bohm-Gross & Landau literature-exact | +| `_spectrax1d` | AW-Hermite $\alpha=\sqrt{2T/m}$ (self-consistent) | **IMMUNE** — zero imports of `normalization.py` | **IMMUNE** | α is a raw config float; physics identical to σ-convention | +| `_hermite_poisson_1d` | AW-Hermite $\alpha=\sqrt{2T/m}$ (self-consistent) | **IMMUNE** — zero imports | **IMMUNE** | Not even registered in ergoExo dispatch | +| `vfp1d` | $\sqrt{2T/m}$ (most-probable, self-consistent) | **IMMUNE** — uses `laser_normalization` ($v_0=c$), never `electron_debye_normalization` | **IMMUNE if** `vth_norm()` is left alone; **BROKEN if** `vth_norm()` is "fixed" | The reason the fix must not touch `vth_norm()` | +| `osiris` wrapper | OSIRIS native ($u_{th}=\sqrt{T/m}/c$, in-deck) | **IMMUNE** — `skin_depth_normalization` ($T_0$=None, $v_0=c$); no eV→uth conversion exists in the wrapper | **IMMUNE** | `vth_norm()` unreachable (would TypeError on $T_0$=None) | +| `vlasov1d2v` | $\sqrt{T/m}$ (σ), $m$≡1 | **IMMUNE by orphan status** — not in ergoExo dispatch, would KeyError before running | **IMMUNE** | Legacy code; carries many un-fixed bugs (§3.9) | + +**Fully immune to both bug and fix:** `_lpse2d`, `_spectrax1d`, `_hermite_poisson_1d`, `osiris`, `vfp1d` (conditional on the fix not touching `vth_norm()`), and `vlasov1d2v` (by virtue of being unrunnable). +**Affected now:** `_vlasov1d`, `_vlasov2d`, `_pic1d` (via `normalization.py`) and `_tf1d` (via its private duplicate). +**Every solver's dimensionless dynamics are unaffected in shipped configs/tests** — the corruption is confined to logged units, physical interpretation of results, dimensional-string inputs, and the EM-branch $\hat c$. + +### 1.2 Why no test ever caught the √2 + +The bug survived every test suite in the repository, and understanding why is itself a finding — the test suites, by construction, cannot see this class of error: + +1. **Every physics test works entirely in dimensionless code units.** Landau damping, ion-acoustic, Bohm-Gross, two-stream, gyro — all set $k$, $\omega$, $T_0$, box sizes as numeric floats, which pass through `normalize()` untouched (`normalization.py:57–58`). The buggy conversion layer is simply never on the tested code path. A purely numeric run *is* a correct simulation under the σ-convention; only its physical labels are wrong, and no test reads the labels. +2. **The theory references the tests compare against live in the same code-unit world.** The kinetic roots come from `electrostatic.py` with `maxwellian_convention_factor=2` and $v_{th}=1$ — i.e., the reference is expressed in the engine's internal convention, so agreement validates internal consistency, not the dimensional dictionary. The tests pin *which* convention the engine uses (invaluable for this audit) but cannot detect that `normalization.py` speaks a different one. +3. **Where a test does touch the dimensional layer, it is self-referential.** `_vlasov2d`'s `test_em_dispersion` computes $\hat c$ from the same `electron_debye_normalization` the solver uses and places the driver at $\omega^2 = 1+\hat c^2 k^2$ with that same $\hat c$ — bug and reference move in lockstep, so it passes with the wrong $c$. +4. **The only tests that pin dimensional values pin the *wrong* ones.** The `_vlasov1d` `*_derived_config.yml` regression fixtures were generated under the buggy convention, so they actively *enforce* the √2 values (`c_light: 11.302`, `x0: 12.14 nm`, …) — a fix makes tests fail, not the bug. +5. **No solver has an absolute physical-units validation** (e.g., a Landau rate checked against a rate in Hz for a stated density and temperature, or an SRS resonance checked against a wavelength in nm), and none has a total-energy conservation diagnostic. Either would have caught the mismatch the first time a dimensional input mattered. + +The general lesson: dimensionless-physics tests validate the engine; only a test that crosses the units boundary — dimensional input in, dimensional observable out, compared against an independent physical reference — can validate the normalization layer. The repo currently has zero such tests. Adding one per normalization entry point (a "round-trip units test") is the cheapest structural guard against recurrence, and belongs in the fix PR alongside item 6 of §4. + +### 1.3 The three thermal-velocity conventions in the codebase + +All are *internally* self-consistent within their modules; the mixing hazard is at module boundaries and in shared helpers: + +| Convention | Definition | Used by | +|---|---|---| +| σ (RMS / std-dev) | $v_{th} = \sqrt{T/m}$, Maxwellian $e^{-v^2/2}$ at $T{=}1$, $L_0=\lambda_{De}$ | `_vlasov1d`, `_vlasov2d`, `_pic1d`, `_tf1d`, `_lpse2d`, `vlasov1d2v` dynamics; `electrostatic.py` (mcf=2 default); OSIRIS `uth` | +| Most-probable | $v_{th} = \sqrt{2T/m}$ | `vfp1d` (with compensating init factor); `vth_norm()`; **`electron_debye_normalization` (the outlier)**; `_tf1d/modules.py:62` (private duplicate, diagnostics-only) | +| AW-Hermite scale | $\alpha = \sqrt{2T/m}$ as basis parameter; represented Maxwellian still has variance $T/m$ | `_spectrax1d`, `_hermite_poisson_1d` (Schumer–Holloway / Parker–Dellar standard) | + +--- + +## 2. Normalization wiring map + +Grep-verified callers of each `normalization.py` entry point: + +| Entry point | Definition | Callers | +|---|---|---| +| `electron_debye_normalization` | $v_0=\sqrt{2T_0/m_e}$ ← **the bug**, $L_0=v_0/\omega_{p0}$ | `_vlasov1d`, `_vlasov2d` (`modules.py:55`), `_pic1d` (`simulation.py:78`) | +| `laser_normalization` | $v_0=c$, $L_0=c/\omega_L$, $n_0=n_{crit}$, $T_0$ not self-consistent with $v_0$ (documented) | `vfp1d` (`base.py:17`) only | +| `skin_depth_normalization(_from_frequency)` | $v_0=c$, $L_0=c/\omega_{p0}$, $T_0$=None | `osiris` only | +| `vth_norm()` | $\sqrt{2T_0/m_0}/v_0$ — hard-codes most-probable | **`vfp1d` only** (grid.py:118, base.py:93, 98, 152) | +| `speed_of_light_norm()` | $c/v_0$ | `_vlasov1d`, `_vlasov2d`, `_pic1d` (√2-tainted); `osiris` (=1, exact) | +| `normalize(s, norm, dim)` | numeric passthrough; strings via $L_0$/τ/$v_0$/$T_0$ | vlasov family, `functions.py`, `vfp1d/helpers.py` profiles | +| — (no normalization import) | | `_tf1d` (own inline pint), `_lpse2d` (own unit system), `_spectrax1d`, `_hermite_poisson_1d` (raw α floats), `vlasov1d2v` (reads cfg keys nothing populates) | + +Shared physics kernels: +- `electrostatic.py` — dispersion/Z-function utilities. Self-consistent under σ-convention: `maxwellian_convention_factor=2` (default, used by **every** caller in the repo) means $f_0\propto e^{-v^2/2v_{th}^2}$, $\xi=\omega/(\sqrt2 k v_{th})$. This is the reference that pins the working convention of `_vlasov1d`, `_vlasov2d`, `_pic1d`, `_tf1d`, `_spectrax1d`, `_hermite_poisson_1d` tests. +- `driftdiffusion.py` — Dougherty/LB kernel; measures $T=\langle(v-\bar v)^2\rangle$ (or spherical $\langle v^4\rangle/3\langle v^2\rangle$), relaxes to that width. Convention-agnostic and self-adapting; verified clean in both audits. +- `functions.py` — envelope/profile layer for the vlasov family. Numeric inputs pass through; dimensional strings go through $L_0$ (√2-tainted until the fix). One bug found (§3.10). +- `utils.py` — no physics content. + +--- + +## 3. Per-module findings + +Finding IDs continue the companion doc's F-series where the bug is a copy; new IDs are per-module. + +### 3.1 `_vlasov2d` — affected now (units/EM), fix-safe, partially cleaned-up lineage + +Convention: σ, identical to `_vlasov1d` (`helpers.py:62` `v_th=√(T0/mass)`, same super-Gaussian α; 2-D init correctly normalized, $\iint f\,d^2v = n$, equal widths in $v_x,v_y$). Landau test pins σ=1 via the shared kinetic roots. Pushers all unity-coefficient; TE-mode Maxwell curl signs verified; initial 2-D Poisson correct; magnetic rotation $\theta=-(q/m)B_z dt$ verified by the gyro test. + +Bug-copy status vs `_vlasov1d`: **F5 fixed** (vbar divides by n, `fokker_planck.py:32`), **F6 fixed** (moments centered on $u=j/n$, `storage.py:100`), **F11 fixed** (`dx=(xmax-xmin)/nx`, `grid.py:60`), F7 N/A (no entropy diagnostic). **F4 still present** (Krook target $e^{-v_x^2/2}e^{-v_y^2/2}$ hard-coded, `fokker_planck.py:127–131`; latent, off in configs). **F10 still present** (`float(cfg.T0/v0x/v0y)`, `simulation.py:158–163`). **F8 present** (same α, latent at m=2). + +√2 exposure: calls `electron_debye_normalization` (`modules.py:55`) → logged units √2-wrong (F2 analog), physical T 2× label (F1 analog), EM $\hat c = c/v_0$ √2-small (F3 analog). `test_em_dispersion` passes *because it is self-consistent, not correct* — it recomputes the same wrong $\hat c$ for the driver. **No regression fixtures exist**, so the fix requires no fixture regeneration; all three tests pass unchanged (the EM test moves in lockstep with the fix). + +New (minor): `mean_KE` has the ½, field proxies $\langle E^2\rangle,\langle B^2\rangle$ don't (no combined conserved-energy diagnostic); `Txx/Tyy` are pressures ($n\cdot$variance) despite the T-name, and `T` omits the mass factor (fine for electrons); the separable per-axis Dougherty never isotropizes $T_x \leftrightarrow T_y$ (no cross-axis coupling operator exists in this module); dead config knobs `UnitsConfig.laser_wavelength`, `Z`. + +### 3.2 `_pic1d` — dynamics immune, logged units affected, fix-safe + +Convention: σ. The particle loaders (`helpers.py:34` quiet inverse-CDF, `:58` random rejection) use `v_thermal = np.sqrt(T0/mass)` — byte-for-byte the vlasov1d initializer convention. Pinned independently by `test_bohm_gross.py:34` (kinetic dielectric with vth=1) and `test_landau_damping.py:38–41` ($\omega^2=1+3k^2$, σ=1 Landau rate). + +Verified clean: deposit/gather use the identical B-spline kernel (momentum-conserving, self-force-free); uniform loading deposits exactly $n=1$ in expectation ($w=n_0L/N$, partition of unity); Poisson $E_k=-i\rho_k/k$ unity prefactor; KDK leapfrog and Yoshida4 push coefficients exact; ponderomotive $(q/m)^2(-\tfrac12\partial_x a^2)$ matches vlasov1d's verified chain. + +√2 exposure: `write_units` (`modules.py:44–57`) logs √2-wrong `v0, x0, c_light, box_length` (a 2000 eV epw run is physically 4000 eV). Dynamics immune: every shipped config and test input is numeric; loader reads raw `T0`. **No fixtures.** Fix is strictly an improvement — nothing breaks. + +New findings: +- **P1 (confirmed):** energy diagnostics mix extensivity: `mean_KE = 0.5·m·Σw v²` is a box integral with the ½; `mean_e2 = mean(e²)` is a per-cell mean without ½ or dx. `mean_KE + mean_e2` is not conserved; no valid energy monitor exists. +- **P2 (confirmed):** quiet and random loaders truncate velocity differently for drifting species — quiet clips to $[-v_{max}, v_{max}]$ absolute (`helpers.py:38`), random clips to $[v_0-v_{max}, v_0+v_{max}]$ (`helpers.py:65`). Same config, two different realized distributions; quiet wrongly clips the high-v side of a drifting Maxwellian. Latent for shipped cold-beam decks. +- **P3 (confirmed):** `mean_KE/mean_p` are sums but named `mean_*` (naming, cf. F6-class). +- Inherited: F11 (`dx=xmax/nx` while particle wrap and placement correctly use $x_{max}-x_{min}$ — the field grid and particle box disagree for $x_{min}\ne0$), F10, F8. + +### 3.3 `_tf1d` — decoupled from `normalization.py`, but carries a private copy of the √2 + +Convention: σ. Derived from the pushers: linearizing continuity + momentum ($-u\partial_x u - \frac{1}{n}\partial_x(p/m) - \frac{q}{m}E$) + adiabatic energy ($\gamma=3$) + Poisson gives $\omega^2 = 1+3k^2$, exactly what `test_resonance.py:19` asserts; the kinetic branch uses the shared roots (mcf=2). The Poisson sign ($E_k=+i\rho_k/k$, i.e. $\partial_x E=-\rho$) and the momentum sign ($-(q/m)E$) are *both* opposite the textbook and cancel — verified stable by derivation and by the passing test. Landau closure decays the field at exactly $2\,\mathrm{Im}\,\omega$, matching `test_landau_damping.py:67`. + +**T1 (confirmed, the headline):** `_tf1d/modules.py:62` re-implements `v0 = √(2·T0/m_e)` inline in its own `write_units()` (with its own pint registry — no import of `normalization.py`). Every tf1d run logs `v0, x0, c_light, beta, box_length, sim_duration` √2-wrong and implies 2× the stated temperature. Damage confined to logged diagnostics (no dimensional-string inputs exist in tf1d; the EM `WaveSolver` branch is commented out; no fixtures). **Fixing `normalization.py` does not fix this** — `modules.py:62` needs its own `2.0*T0 → T0` in lockstep. + +Other findings: +- **T2 (confirmed):** the kinetic-γ pressure closure (`pushers.py:195` `wr_corr=(wrs²-1)/k²` + γ forced to 1) yields $\omega^2 = (w_{rs}^2-1)T_0+1$ — exact only at $T_0=1$. Convention- and T₀-locked, same class as vlasov1d's Krook (F4). Latent (all configs use $T_0=1$). +- **T3 (confirmed):** `modules.py:78` `nuee_norm = nuee/wp0` — the *correct* Hz→code conversion — is computed and then **discarded** (never logged, never used). Meanwhile `physics..trapping.nuee` is an unrelated raw ML-input float. There is no collisional friction in the tf1d equations at all. (F9-class trap.) +- **T4 (confirmed):** `docs/source/usage/tf1d.md` momentum equation has a spurious $1/n$ on the E-force and mislabels the Poisson equation. Docs-only. +- **T5 (suspected, latent):** `EnergyStepper` evolves $p$ but advects/compresses $p/m$ — ambiguous mass normalization of the pressure variable; inconsistent with the momentum equation for mobile ions ($m=1836$). All configs have ions off. +- **T9 (minor):** `resonance_search.yaml` sets `ion.landau_damping: True` with `ion.is_on: False` (inert, misleading). + +### 3.4 `_lpse2d` — fully immune; thermally clean; unrelated bugs found + +Physical-units solver (ps/µm, Gaussian-cgs-derived) with its own `write_units` (`helpers.py:71`); zero imports of `normalization.py`. Convention: $v_{te} = c\sqrt{T_e/511}$ = $\sqrt{T_e/m_e}$ (σ) used consistently in the Bohm-Gross term ($e^{-i\,1.5\,v_{te}^2 k^2/\omega_{pe}\,dt}$, coefficient 3/2 correct for σ), the Landau damping rate (verified literature-exact: $\sqrt{\pi/8}(1+\tfrac32 k^2\lambda_D^2)(k\lambda_D)^{-3}e^{-3/2-1/(2k^2\lambda_D^2)}$), the sound speed, and the TPD threshold. Numerically cross-checked: the test config's driver $\omega = 1.5k^2v_{te}^2/\omega_{p0} = 19.7 \approx 20$ as configured. TPD/SRS coupling constants contain no thermal factor (convention-independent); laser $E_0=\sqrt{8\pi I/c}$ and WKB swelling $(1-n/n_c)^{-1/4}$ verified. + +Findings (none √2-class): +- **L1 (confirmed):** logged `lambda_D = vte/w0` (`helpers.py:111`) divides by the **laser** frequency instead of $\omega_{pe}$ — factor 2 too small at $n_c/4$. Diagnostic only (dynamics compute $\omega_{pe}^2/k^2v_{te}^2$ directly). +- **L2 (confirmed, latent runtime bug):** the `E2` electrostatic-driver path calls `self.epw.driver(...)` (`core/vector_field.py:84`) but the active `SpectralEPWSolver` has no `driver` attribute (only the commented-out `SpectralPotential` does) → `AttributeError` for any config with an `E2` driver. Latent only because `test_epw_frequency` is currently disabled (`pass` body) — meaning **the module's dispersion conventions have no live regression test**. +- **L3 (confirmed):** `core/trapper.py` is dead code (never instantiated by `SplitStep`); it is a stale duplicate of `_tf1d`'s trapper. Note: it uses `electrostatic.py`, *not* `driftdiffusion.py` — the vlasov1d audit's cross-module note claiming lpse2d uses the Dougherty kernel is wrong for the current tree. +- **L4 (confirmed):** the σ-convention is undocumented in the lpse2d docs — the main *future* √2 risk here is someone "harmonizing" it onto `normalization.py`. +- **L5 (minor):** the Landau/Bohm-Gross formulas are duplicated byte-for-byte in `SpectralPotential` and `SpectralEPWSolver` (drift risk); `nu_coll`'s `/2` (amplitude vs energy rate) is correct but easy to double-count. + +### 3.5 `_spectrax1d` and `_hermite_poisson_1d` — fully immune; self-consistent AW-Hermite basis + +Both use the asymmetrically-weighted Hermite basis (Schumer–Holloway / Parker–Dellar) with scale $\alpha = \sqrt{2T/m}$ taken as a **raw config float** — zero imports of `normalization.py`. The represented Maxwellian has variance $\alpha^2/2 = T/m$: physically identical to the σ-convention solvers. Verified by derivation (streaming ladder $\alpha\sqrt{n/2}$ + force ladder $\sqrt{2n}/\alpha$ + Ampère/Poisson coupling → Langmuir $\omega=\omega_{pe}$ exactly, Bohm-Gross $\omega^2=1+3(k\lambda_{De})^2$ with $\lambda_{De}=\alpha/\sqrt2$) and by tests: both modules' Landau tests set $\alpha_e = \sqrt2\,k\lambda_D/k$ explicitly against the shared mcf=2 kinetic roots (2%/5% tolerance for hermite-poisson; a cross-module test pins the two modules' E-coupling equal to 1e-12). + +The coincidence that $\alpha=\sqrt{2T/m}$ matches the `normalization.py:91` outlier is harmless — there is no code coupling in either direction. Both integrator paths (DoPri8 and Lawson-RK4 exponential operators) share identical ladder coefficients. The historical inverted E-coupling bug in hermite_poisson (`C[n+1]` vs `C[n-1]`) is fixed and regression-locked by three tests. + +Findings (all diagnostic-label/minor): **S1** logged `lambda_D` is the *total* (ion-dominated) Debye length while `k_norm` uses electron-only and hard-codes mode 1 (`base_module.py:186,189`); **S2** the ion "temperature" diagnostic is a velocity variance missing the $m_i$ factor ($T_i/m_i$, misleading next to $T_e$; `storage.py:386–405`); **S3** shipped `landau-damping.yaml` has `Nn: 4` — far too few Hermite modes to resolve the damping it is named for (tests override to 512/32); `_spectrax1d/helpers.py` is an all-stub file not on any code path. + +### 3.6 `vfp1d` — immune to the bug; the reason the fix must not touch `vth_norm()` + +Uses `laser_normalization` (`base.py:17`): $v_0 = c$, $L_0 = c/\omega_L$, $n_0 = n_{crit}$, $T_0$ = reference eV (documented as not self-consistent with $v_0$). Velocity grid in units of $c$. Thermal convention: **most-probable, $v_{th} = \sqrt{2T/m}$**, deliberately and self-consistently: + +- `vth_norm()` supplies $\sqrt{2T_0/m}/c$ (4 call sites, all vfp1d); +- the initializer's width factor `helpers.py:111` is $\alpha = \sqrt{3\Gamma(3/m)/(2\Gamma(5/m))}$ — note the extra `/2` vs the vlasov-family α — giving $\alpha=1$ at $m{=}2$ and realized variance exactly $T/m$ (the √2 and the /2 cancel); +- the v-grid extent `grid.py:118` `vmax = 8·vth_norm()/√2` = 8 standard deviations (the /√2 converts most-probable→σ); +- the IB coefficient `base.py:81` ($0.093373\,\lambda_{\mu m}^2/T_{keV}$ per $10^{15}$ W/cm²) normalizes $v_{osc}^2$ to $2T/m$, consistent; +- the temperature diagnostic `storage.py:317–322` uses the spherical variance $T = \langle v^4\rangle/3\langle v^2\rangle$ — reads the true temperature, no √2. + +Collision operators ($\nu_{ee}$ coefficient anchored to $v_0=c$ via $r_e$, Lorentz e–i with Epperlein–Haines Z*, Rosenbluth I/J integrals) are temperature-convention-independent. **No internal √2 inconsistency exists.** + +**C1 (confirmed, cross-module, high):** because items above are keyed to `vth_norm()` while the compensating α is a separate constant, redefining `vth_norm()` → $\sqrt{T_0/m}/v_0$ (as the vlasov1d audit's §5.1 originally suggested) silently initializes VFP-1D at **half** the intended temperature while all transport coefficients stay at $T_0$; the Spitzer/Epperlein–Haines gold test (`test_kappa_eh.py`, $\kappa\propto T^{5/2}$) would likely fail, and the pure-operator tests would *not* catch it. **The fix must be confined to `electron_debye_normalization:91–92`.** (The companion doc's §5.1 has been amended accordingly.) + +Other findings: **C2** `base.py:118` stores the bound method `norm.vth_norm` instead of calling it (harmless repr-string in cfg; wrong); **C3** `storage.py:170` hard-codes $9.09\times10^{21}$ cm⁻³ per $n_c$ — that is $n_{crit}(351\,\mathrm{nm})$; wrong labeling for any other `laser_wavelength` (should derive from `norm.n0`); **C4** = F11 copy (`grid.py:71` `dx=xmax/nx`); **C6/C7 (suspected, latent — IB not in any shipped config):** the production inverse-bremsstrahlung wiring looks under-normalized — `w0_norm` is identically 1.0 under laser normalization, the Langdon-factor argument $Z^2 n_i/(\omega_0 v^3)$ carries no collision-frequency coefficient, and `vosc2_per_intensity` (thermal-speed² units) is consumed as $v_{osc}^2$ in grid ($c$) units — a $\sim(c/v_{th})^2$ discrepancy between the unit tests (where the grid unit *is* the thermal speed) and production. Recommend a Spitzer-IB validation run before enabling IB in production. + +### 3.7 `osiris` wrapper — fully immune + +Uses `skin_depth_normalization(_from_frequency)` ($v_0=c$, $T_0$=None). The wrapper drives a native OSIRIS deck: **no eV→`uth` conversion exists anywhere in the wrapper** — `uth` values pass through verbatim from the deck, so the T→v √2 trap cannot occur. Verified exact: $a_0 \to E_{peak} = a_0\omega_L m_e c/e \to I = E^2\epsilon_0 c/2$ (matches $I\lambda_{\mu m}^2 = 1.37\times10^{18}a_0^2$, test-confirmed); `c_light = beta = 1.0` identically; box lengths in true skin depths; `units.yaml` deliberately omits temperature-dependent keys (consistent with $T_0$=None). Temperature diagnostics (`plots.py:1004` $T=\sum u_{th,i}^2$ in $m_ec^2$ units; `:1059` momentum variance) carry no spurious factors. `vth_norm()` is unreachable (would TypeError on $T_0$=None). + +**C8 (low, latent):** the adaptive-box feature's `reference_density = 0.25` "quarter-critical" (`density.py:75`) assumes the deck's $n_0 = n_{crit}$ (true for LPI decks with `omega0=1`); a deck normalized to a different density would silently mis-scale the density bounds (the *length* normalization stays correct). + +### 3.8 Shared `electrostatic.py` and `functions.py` + +`electrostatic.py`: self-consistent under the σ-convention throughout. $Z$ via Faddeeva, $Z' = -2(1+\xi Z)$ correct; `maxwellian_convention_factor=2` (default, used by every caller in the repo) ⇒ $f_0\propto e^{-v^2/2v_{th}^2}$, root returned as $\xi k v_{th}\sqrt{2} = \omega$. One cosmetic blemish (**T7**): the Newton `initial_root_guess` is ω-scaled ($\sqrt{\omega_p^2+3k^2v_{th}^2}$) but the root variable is ξ-scaled — converges anyway. + +`functions.py`: numeric inputs pass through; `EnvelopeFunction` amplitudes (`baseline`, `bump_height`) are raw floats (the F9 collision-frequency trap); `SineFunction.wavenumber` correctly uses `dim="k"`. **T8 (suspected):** `LinearFunction`/`ExponentialFunction` normalize `val_at_center` — a *density* — with `dim="x"` (`functions.py:159–161, 178–180`); dimensionally wrong for any string input. + +### 3.9 `vlasov1d2v` — orphaned legacy; a museum of the known bugs plus three new ones + +Not registered in the ergoExo dispatch (`_base_.py:294–333` → NotImplementedError); reads `cfg["units"]["derived"]` keys nothing populates → cannot run. Never calls `normalization.py`. Convention: σ (init `exp(-(v_x^2+v_y^2)/2T_0)`, mass≡1; the super-Gaussian machinery is commented out so `m` is silently ignored). + +Carries live copies of: **F5** (`fokker_planck.py:93,98` — vbar *and* the diffusion coefficient computed without ÷n), **F6 + F7** (`storage.py:152–158, 296`), **F10**, **F11** (`helpers.py:211`). F4 present but dead (Krook never instantiated). + +New bugs (would matter if ever revived): +- **N1 (confirmed):** the density profile is computed for every basis and folded into the ion background (`helpers.py:287`), but the line applying it to $f$ (`helpers.py:83–85`) is **commented out** — electrons always start uniform, so non-uniform configs begin with spurious net charge and no actual density perturbation in $f$. +- **N2 (confirmed):** the explicit Dougherty substeps apply ν twice (`fokker_planck.py:134–135, 141–142`: `dfdt = nu*ddx(...)` then `f + dt*nu*dfdt`) — collisional relaxation runs at $\nu^2$. +- **N4 (confirmed):** `integrator.py:234` passes the **nu_ee** envelope args when building the **nu_ei** profile (copy-paste) — the configured nu_ei profile is silently ignored. +- **N3 (design note):** without collisions the $v_y$ dimension is dynamically inert (no $v_y$ force, streaming by $v_x$ only) — the entire $n_{vy}$ grid is dead compute except through the FP operators. +- Minor: default scalars broadcast over the $v_y$ axis only (`storage.py:292–295`). + +Recommendation: either delete `vlasov1d2v` or quarantine it with a module-level comment; in its current state it is a trap for anyone who greps for reference implementations. + +### 3.10 Cross-module systemic observations + +1. **Driver amplitude conventions differ across the Vlasov family** (N5): `_vlasov1d` Ex driver $\propto \omega a_0$; `_vlasov2d` current source $\propto \omega^2 a_0$; `vlasov1d2v` $\propto |k| a_0$. Same config key, three physical meanings. Worth unifying or documenting per-solver. +2. **No solver has a valid total-energy conservation diagnostic.** vlasov1d/vlasov2d/pic1d all mix ½-factors and extensivity between kinetic and field terms. A per-solver conserved-energy scalar is the single cheapest guard that would have caught the √2 class of bugs. +3. **The F-series bugs are lineage-correlated:** `_vlasov2d` (newest refactor) fixed F5/F6/F11; `_pic1d` inherits F10/F11 by importing `_vlasov1d`'s grid/simulation; `vlasov1d2v` (oldest) has everything. When fixing `_vlasov1d`, fix the copies in the same PR: Krook targets (`_vlasov1d`, `_vlasov2d`, dead `vlasov1d2v`), `dx`/xmin (`_vlasov1d/grid.py:55`, `vfp1d/grid.py:71`, `vlasov1d2v/helpers.py:211`), raw-float T0/v0 (all three vlasov-family simulation.py files). +4. **Hz→code-unit collision-frequency conversion** exists correctly in exactly one place (`_tf1d/modules.py:78`) and is dead there; nowhere is it live. Every solver that takes a collision frequency takes it as a raw code-unit float with undocumented units. + +--- + +## 4. Consolidated fix checklist (supersedes §5 of the vlasov1d audit where they differ) + +**The √2 fix, safely scoped:** +1. `normalization.py:91–92`: $v_0 \to \sqrt{T_0/m_e}$, $x_0 \to \lambda_{De}$ — inside `electron_debye_normalization` **only**. +2. **Do NOT change `vth_norm()`** (`normalization.py:48–50`) — it is vfp1d-private and correct for vfp1d (C1). If its name is judged misleading, rename to `most_probable_speed_norm()`; do not change the value. +3. `_tf1d/modules.py:62`: `2.0*T0 → T0` (private duplicate, T1). +4. Regenerate the three `_vlasov1d` `*_derived_config.yml` regression fixtures (only vlasov1d has fixtures). +5. Fix the convention-locked Krook targets in the same PR (`_vlasov1d/…/fokker_planck.py:245`, `_vlasov2d/…/fokker_planck.py:127–131`): build from species T0/mass. +6. Re-validate: vlasov1d Landau + ion-acoustic (must pass unchanged), vlasov2d Landau/gyro/EM-dispersion (EM moves in lockstep), pic1d all three (unchanged), vfp1d `test_kappa_eh` (must pass unchanged — the sentinel that `vth_norm()` was left alone), numeric-driver SRS configs re-tuned per vlasov1d F3. +7. Document the per-solver conventions in each solver's docs page (σ vs most-probable vs AW-α; see §1.3) — the biggest residual risk is a future "harmonization" that imports the wrong convention into a currently-clean module (especially `_lpse2d`, L4). + +**Independent bug fixes, by priority:** +- High (live, wrong physics if exercised): vfp1d IB normalization audit before production use (C6/C7); lpse2d `E2` driver path + re-enable `test_epw_frequency` (L2). +- Medium (live, wrong diagnostics): lpse2d logged `lambda_D` (L1); pic1d/vlasov2d/vlasov1d energy-diagnostic normalization + add conserved-energy scalars; spectrax ion-T mass factor (S2); vfp1d hard-coded $n_{crit}(351\,nm)$ (C3); pic1d quiet-loader drift truncation (P2). +- Low / hygiene: tf1d docs (T4) and dead `nuee_norm` (T3 — log it instead of dropping it); `functions.py` `val_at_center` dim (T8); vfp1d bound-method store (C2); spectrax `Nn:4` config (S3); delete or quarantine `vlasov1d2v` and lpse2d's dead trapper (L3); unify driver-amplitude conventions or document them (N5). + +--- + +## 5. What was verified clean (highlights) + +- All solver cores' dimensionless coefficients: vlasov1d/2d pushers and field solves, pic1d push/deposit/gather/Poisson, tf1d fluid system (including its double sign cancellation, verified by derivation), Hermite ladder algebra in both spectral modules (verified by derivation to the Langmuir and Bohm-Gross limits), lpse2d Bohm-Gross/Landau/TPD/SRS coefficients against literature, vfp1d collision coefficients and temperature diagnostics, osiris a0↔intensity conversions. +- The shared `driftdiffusion.py` Dougherty kernel (both audits): self-adapting, Buet factor-2 correctly absorbed, spherical temperature correct. +- `electrostatic.py`'s Z-function algebra and mcf=2 semantics — the single reference that consistently pins six modules' tests to the σ-convention. + +--- + +*Companion deep-dive for `_vlasov1d`: `VLASOV1D_CONVENTIONS_AUDIT.md`. Method: six parallel independent module audits (one Opus agent per slice) with all load-bearing claims re-verified against source before inclusion. No source files were modified.* diff --git a/docs/audits/2026-07-16_VLASOV1D_CONVENTIONS_AUDIT.md b/docs/audits/2026-07-16_VLASOV1D_CONVENTIONS_AUDIT.md new file mode 100644 index 00000000..81a302a5 --- /dev/null +++ b/docs/audits/2026-07-16_VLASOV1D_CONVENTIONS_AUDIT.md @@ -0,0 +1,249 @@ +# Vlasov-1D Physics & Conventions Audit + +**Date:** 2026-07-16 +**Scope:** `adept/_vlasov1d/` (all files), plus the shared code it depends on: `adept/normalization.py`, `adept/driftdiffusion.py`, `adept/electrostatic.py`, and the `tests/test_vlasov1d/` suite and example configs. +**Trigger:** the discovery that `normalization.py:91` uses the $v_0 = \sqrt{2T_0/m}$ (most-probable-speed) convention while `_vlasov1d/helpers.py:69` uses $v_{th} = \sqrt{T_0/m}$ (RMS / standard-deviation convention). Neither has been changed yet; this audit determines which convention the module actually runs in, inventories every equation, and lists all inconsistencies found. + +--- + +## 1. Executive summary + +**The engine's working convention is $v_0 = \sqrt{T_0/m}$.** The distribution initializer, the Fokker–Planck operator, the Krook operator, the Landau-damping and ion-acoustic tests (the module's quantitative gold standards), and every numeric example config are all mutually consistent under: code velocity in units of the thermal *standard deviation* $\sqrt{T_0/m_e}$, code length unit $L_0 = \lambda_{De}$, code wavenumber $= k\lambda_{De}$, Maxwellian $\propto e^{-v^2/2}$ at $T=1$. + +**`adept/normalization.py:91–92` is the sole outlier** (`electron_debye_normalization`: $v_0 = \sqrt{2T_0/m_e}$, $L_0 = \sqrt{2}\,\lambda_{De}$). Because numeric config inputs pass through `normalize()` untouched, this does **not** corrupt the dynamics of numeric-input runs — but it means: + +- A run with `normalizing_temperature: 2000eV` and species `T0: 1.0` is physically a **4000 eV** plasma (factor 2 in temperature). +- Every dimensional **string** input converted with $L_0$ (box sizes in µm, gradient scale lengths, laser $k_0$) is off by $\sqrt{2}$. +- Every **logged** physical unit (`v0`, `x0`, `c_light`, `box_length`) is off by $\sqrt{2}$, and the regression fixtures currently lock in those wrong values. +- The EM wave speed $\hat c = c/v_0$ fed to the wave solver is $\sqrt2$ too small relative to the engine's thermal unit (with a partial cancellation for wavelength-specified drivers; see F3). + +**Recommended fix direction:** change `normalization.py` to $v_0 = \sqrt{T_0/m_e}$ (so $L_0 = \lambda_{De}$), *not* the initializer. Fixing `helpers.py:69` instead (to $\sqrt{T_0/2m}$) would break the currently-passing Landau-damping and ion-acoustic tests and invalidate every existing config's driver $(k_0, \omega_0)$ values. Section 5 lists everything that must move in lockstep. + +The collisionless solver core (`vector_field.py`, pushers) was verified **clean**: every coefficient is exactly unity under the declared normalization and is convention-independent. The convention bug lives entirely in the dimensional-conversion layer. + +Beyond the $\sqrt2$ issue, the audit found several independent bugs and traps, listed in Section 4 (notably: FP `compute_vbar` missing $1/n$; storage central moments centered on $n u$ instead of $u$; a sign inconsistency between the two `-flogf` entropy diagnostics; the super-Gaussian $\alpha$ fixing the wrong moment for $m\ne2$; the Krook target hard-coded to $T=1$; and the collision-frequency normalization chain being entirely manual). + +--- + +## 2. The two candidate conventions + +With $\omega_{p0} = \sqrt{n_0 e^2/(\epsilon_0 m_e)}$, $\tau = 1/\omega_{p0}$ in both cases: + +| Quantity | (a) `normalization.py` as written | (b) engine's actual convention | +|---|---|---| +| Velocity unit $v_0$ | $\sqrt{2T_0/m_e}$ (most-probable speed) | $\sqrt{T_0/m_e}$ (RMS / std-dev) | +| Length unit $L_0 = v_0/\omega_{p0}$ | $\sqrt{2}\,\lambda_{De}$ | $\lambda_{De}$ | +| Maxwellian at $\hat T = 1$ | $\propto e^{-v^2}$ (variance ½) | $\propto e^{-v^2/2}$ (variance 1) | +| Code wavenumber $\hat k$ | $\sqrt2\, k\lambda_{De}$ | $k\lambda_{De}$ | +| Bohm–Gross | $\hat\omega^2 = 1 + \tfrac{3}{2}\hat k^2$ | $\hat\omega^2 = 1 + 3\hat k^2$ | +| $\hat c = c/v_0$ (2000 eV) | 11.30 | 15.98 | + +**Why the test suite never caught this:** every physics test (Landau, ion-acoustic, EM dispersion) works in dimensionless code units — numeric inputs bypass `normalize()` entirely (`normalization.py:57–58`), so the buggy conversion layer is never on the tested path, and the theory references (`electrostatic.py`, mcf=2, $v_{th}=1$) are expressed in the engine's internal convention, so agreement validates internal consistency rather than the physical-units dictionary. Worse, the only tests that *do* pin dimensional values — the `*_derived_config.yml` regression fixtures — were generated under the buggy convention and therefore lock the wrong values in. No test crosses the units boundary (dimensional input in, dimensional observable out, against an independent physical reference), which is precisely the kind of test that must accompany the fix (see `ADEPT_CONVENTIONS_AUDIT.md` §1.2 for the full analysis). + +Evidence pinning convention (b) as the working one, strongest first: + +1. **Landau-damping test** (`tests/test_vlasov1d/test_landau_damping.py:24`) calls `electrostatic.get_roots_to_electrostatic_dispersion(1.0, 1.0, k)` — i.e. $\omega_{pe}=1$, $v_{th}=1$ with `maxwellian_convention_factor=2` (`adept/electrostatic.py:79,115,98`), which is the kinetic dielectric for $f_0 \propto e^{-v^2/2v_{th}^2}$, $\xi = \omega/(\sqrt2 k v_{th})$. The measured $\omega(k)$ matches these roots to 2 decimals only if the simulated Maxwellian has $\sigma = 1$. (E.g. $k=0.3 \Rightarrow \omega = 1.1598$, which is exactly the `w0` in `resonance.yaml` and `epw.yaml`; under convention (a) the resonance would sit near 1.07 and the test would fail.) +2. **Ion-acoustic test** (`test_ion_acoustic_wave.py`) uses $c_s^2 = ZT_e/m_i$ and $\omega^2 = k^2c_s^2/(1+k^2\lambda_D^2)$ with $\lambda_D = 1$ in code units — asserting $L_0 = \lambda_{De}$. +3. **Initializer** (`helpers.py:69–76`), **Krook target** (`fokker_planck.py:245`), and **Dougherty stationary state** (`fokker_planck.py:51,96` + `driftdiffusion.py`) all build/assume $e^{-v^2/(2T/m)}$, variance $T/m$. +4. Every numeric example config (`epw.yaml`, `resonance.yaml`, `wavepacket.yaml`, `bump-on-tail.yaml`) uses driver $(k_0,\omega_0)$ pairs consistent only with $\hat k = k\lambda_{De}$ and $\sigma = \sqrt{T/m}$. + +--- + +## 3. Equation inventory + +### 3.1 Normalization layer (`adept/normalization.py`) + +| Location | Expression | Meaning | Convention note | +|---|---|---|---| +| :88 | $\omega_{p0} = \sqrt{n_0e^2/\epsilon_0 m_e}$, $\tau = 1/\omega_{p0}$ | time unit | standard | +| **:91** | $v_0 = \sqrt{2T_0/m_e}$ | velocity unit | **outlier — convention (a)** | +| :92 | $x_0 = v_0/\omega_{p0} = \sqrt2\lambda_{De}$ | length unit | inherits $\sqrt2$ | +| :48–50 | `vth_norm()` $= \sqrt{2T_0/m_0}/v_0 = 1$ | "thermal velocity" | convention (a); **never called** in `_vlasov1d` | +| :52–53 | `speed_of_light_norm()` $= c/v_0$ | $\hat c$ | inherits $\sqrt2$ (F3) | +| :37–38 | NRL $\log\Lambda_{ee}$ | Coulomb log | uses reference $T_0$ | +| :45 | $\nu_{ee} = 2.91\times10^{-6}\, n_{cc}\log\Lambda\, T_{eV}^{-3/2}$ Hz | NRL e–e rate | **diagnostic only**, never converted to code units (F9) | +| :56–74 | `normalize()`: x/L0, t/τ, v/v0, T/T0, k·L0 | dim→code | numeric inputs pass through untouched; `dim="temp"` and `dim="v"` branches are dead code for `_vlasov1d` | + +### 3.2 Initialization & grids (`helpers.py`, `grid.py`, `modules.py`, `simulation.py`, `datamodel.py`) + +| Location | Expression | Meaning | Convention note | +|---|---|---|---| +| helpers.py:65–66 | `dv = 2 vmax/nv`; cell-centered `vax` on $[-v_{max}+dv/2,\ v_{max}-dv/2]$ | v-grid | edges at $\pm v_{max}$; half-cell center offset | +| **helpers.py:69** | $v_{th} = \sqrt{T_0/m}$ | thermal width | **convention (b)** | +| helpers.py:72 | $\alpha = \sqrt{3\,\Gamma(3/m)/\Gamma(5/m)}$ ($=\sqrt2$ at $m{=}2$) | super-Gaussian width factor | fixes $\langle v^4\rangle/\langle v^2\rangle = 3v_{th}^2$; variance $=v_{th}^2$ **only** at $m=2$ (F8) | +| helpers.py:74–76 | $f \propto \exp[-|(v-v_d)/(\alpha v_{th})|^m]$ | init EDF; $m{=}2$: $e^{-(v-v_d)^2/(2T_0/m)}$ | variance $T_0/m$ — convention (b) | +| helpers.py:80,84 | `f /= sum(f)·dv`; `f *= n_prof` | $\int f\,dv = n(x)$ | midpoint rule, consistent with all moments | +| grid.py:55 | `dx = xmax/nx` | cell width | **ignores `xmin`** (F11) | +| grid.py:59–60 | `dt = min(dt, 0.95·dx/c)`, $c = 1/\beta$ | EM CFL | $\hat c$ inherits F1's $\sqrt2$ | +| grid.py:65–66,75 | `nt = int(tmax/dt+1)`; `t = linspace(0, dt·nt, nt)` | time axis | overshoot ≤ dt; `grid.t` spacing ≠ dt (F12) | +| grid.py:74,77 | cell-centered x; `kx = 2π·fftfreq(nx, d=dx)` | spectral grid | inherits F11's dx for `xmin≠0` | +| modules.py:112 | `box_length = (xmax−xmin)·L0` | box size in µm | uses `xmin` correctly (unlike grid.py:55); value carries F1's $\sqrt2$ | +| simulation.py:213–214 | `v0 = float(cfg.v0)`; `T0 = float(cfg.T0)` | species drift & temperature | bare code floats, bypass `normalize()` (F10) | +| simulation.py:74–76 | $a_0 = a_{0,std}\cdot(c/v_0)$ | intensity → quiver velocity in $v_0$ units | internally consistent with pusher (§3.3) | +| simulation.py:79–85 | $\hat k_0 = k_{phys}L_0$, $\hat\omega_0 = \omega_{phys}\tau$ | wavelength driver → code units | $\hat k_0$ carries F1's $\sqrt2$; $\hat\omega_0/\hat k_0 = \hat c$ ✓ | +| simulation.py:87 | `dw0 = 0.0 # ???` | frequency offset placeholder | unresolved (F12) | + +### 3.3 Solver core (`solvers/vector_field.py`, `solvers/pushers/vlasov.py`, `solvers/pushers/field.py`) — **verified clean** + +The dimensionless system implemented, with every coefficient exactly unity under the declared normalization (and independent of the (a)/(b) choice, since $v_0$ cancels): + +$$\partial_t f_s + v\,\partial_x f_s + \frac{\hat q_s}{\hat m_s}\left(E - \frac{\hat q_s}{2 \hat m_s}\partial_x a^2\right)\partial_v f_s = C[f_s]$$ +$$\partial_x E = \sum_s \hat q_s \int f_s\,dv \ (+\ \text{static ion background}), \qquad \partial_t E = -\sum_s \hat q_s \int v f_s\,dv$$ +$$\partial_t^2 a = \hat c^2\,\partial_x^2 a - n_e\, a + S$$ + +| Location | Code | Check | +|---|---|---| +| vlasov.py:210–214, 220 | exact spectral shift $f(x - v\,dt)$ | coefficient $v_0\tau/L_0 = 1$ ✓ | +| vlasov.py:83–87, 127–129 | `force = q·e + (q²/m)·pond; accel = force/m` | $\hat q/\hat m$ E-push; ponderomotive $(\hat q^2/\hat m^2)$ ✓ | +| vector_field.py:411 | `pond = −0.5·∂ₓ(a²)` | exact instantaneous $-\frac{1}{2m}\partial_x p_\perp^2$ with $p_\perp = qA$; the $c/v_0$ scaling of $a_0$ (simulation.py:76) makes $a$ the quiver velocity in $v_0$ units — **no hidden $\sqrt2$; verified exact** | +| field.py:203–224 | $E_k = -i\rho_k/k$ | Poisson prefactor $e^2n_0/(\epsilon_0 m_e\omega_{p0}^2) = 1$ ✓ | +| field.py:263–264, 280 | $E^{n+1} = E^n - dt\,j$ | Ampère prefactor 1 ✓ | +| field.py:337–343 | $\Delta E_k = -i\frac{q}{k}\int dv\, f_k(e^{-ikv\,dt}-1)$ | exact (Hamiltonian) Ampère along free streaming ✓ | +| field.py:146–153 | leapfrog wave eq., plasma term $-n_e a$ | dispersion $\hat\omega^2 = \hat n_e + \hat c^2\hat k^2$ ✓; $n_e$ is electron-only (~$m_e/m_i$ approx., F12) | +| vector_field.py:292–301, 341 | $n_e = -\,\hat q_e \int f_e\,dv$, time-centered | sign keyed to electron charge $-1$ ✓ | +| field.py:21–26 | $E_x$ driver $= (\omega_0{+}\delta\omega)\,a_0\sin(k_0x - \omega t)$ | vector-potential amplitude convention ($E = -\partial_t A$); matches docs (config.md:312) ✓ | +| field.py:67–69, 78–80 | point source $F_0 = 2\omega\hat c\,a_0$ | reproduces amplitude $a_0$ in vacuum; plasma value larger by $k_{vac}/k_{plasma}$, documented in docstring ✓ | + +Charge-sign conventions were checked across Poisson / Ampère / wave-equation and the force term: consistent, no compensating double-error. Current density correctly has **no** mass factor in code (but see docstring bug F12). + +### 3.4 Collisions (`solvers/pushers/fokker_planck.py`, `adept/driftdiffusion.py`) + +Operator: $\partial_t f = \nu\,\partial_v[(v - \bar v)f + T\,\partial_v f]$ (Lenard–Bernstein/Dougherty, Buet notation $\beta = 1/2T$, $D = 1/2\beta = T$, drift $C = 2\beta D(v-\bar v) = (v-\bar v)$). + +| Location | Code | Check | +|---|---|---| +| driftdiffusion.py:101ff (`discrete_temperature`) | $T = \sum f(v-\bar v)^2 dv \,/\, \sum f\,dv$ | full 2nd central moment, no factor 2, no $m$ — variance convention, matches init ✓ | +| driftdiffusion.py:170,180 | $\beta_{init} = 1/2T$; $f_{mx} = e^{-\beta(v-\bar v)^2}$ | stationary state $e^{-(v-\bar v)^2/2T}$ ✓ | +| driftdiffusion.py:333–334 | $C = 2\beta D\,(v_{edge} - \bar v)$ | $2\beta D = 1$ exactly ✓ | +| driftdiffusion.py (flux/Chang–Cooper, implicit solve) | conservative central & positivity-preserving fluxes; $(I - dt\,\nu L)f^{n+1} = f^n$ | numerics only; $\nu\,dt$ dimensionless ✓ | +| driftdiffusion.py:20–25 | Buet "extra factor of 2" note | correctly absorbed into $\beta = 1/2T$; **no stray factor survives** (verified) | +| **fokker_planck.py:81** | `compute_vbar` $= \sum f\,v\,dv$ | returns $n\bar u$, **not** $\bar u$ — missing $1/n$ (F5) | +| **fokker_planck.py:245** | Krook target $f_{mx} \propto e^{-v^2/2}$ | hard-coded variance 1 = ($T{=}1$, $m{=}1$, convention (b)) — convention-locked (F4) | +| fokker_planck.py:262–266 | $f \to f e^{-\nu dt} + n(x) f_{mx}(1 - e^{-\nu dt})$ | BGK; conserves $n$, not momentum/energy (by design) | + +Key property: the **Dougherty operator is self-adapting** — it measures $T = \langle(v-\bar v)^2\rangle$ from $f$ and relaxes toward exactly that width, so it preserves whatever convention the initializer used and needs no change under a convention fix. The **Krook operator does not** — its width is a compile-time constant (F4). + +### 3.5 Diagnostics & storage (`storage.py`) + +| Location | Code | Meaning | Note | +|---|---|---|---| +| :134–135 | $\int(\cdot)\,dv$ = `sum·dv` | midpoint rule | consistent with init normalization ✓ | +| :139 | `n` $= \int f\,dv$ | density | ✓ | +| :140 | `v` $= \int v f\,dv$ | **raw first moment $= n\bar u$**, labeled "v" | flux, not velocity (F6) | +| :141–142 | `p` $= \int (v - n\bar u)^2 f\,dv$ | "pressure" | centered on $n\bar u$, not $\bar u$ (F6); for $n{=}1$ Maxwellian at $T_0$: `p` $= T_0$ under convention (b) ✓, would read $2T_0$ under (a) | +| :143 | `q` $= \int (v-n\bar u)^3 f\,dv$ | heat flux | same centering issue | +| :144 | `-flogf` $= \int f\log|f|\,dv$ | "entropy" | **missing minus sign** (F7) | +| :307 | `mean_-flogf` $= \langle\int -|f|\log|f|\,dv\rangle$ | entropy | has the minus; opposite sign to :144 (F7) | +| :303–306, 308 | `mean_P` $=\langle\int v^2f\rangle$, `mean_j` $=\langle\int vf\rangle$, `mean_n`, `mean_q` $=\langle\int v^3 f\rangle$, `mean_f2` | raw moments | `mean_j` and field "v" are the same integrand under different names (F6) | +| :311–312 | `mean_de2` $=\langle de^2\rangle$, `mean_e2` $=\langle e^2\rangle$ | field energy proxies | no ½; kinetic `mean_P` also lacks ½, consistently — internal conservation OK, but no combined energy monitor exists (F12) | +| :154, :313 | `pond` $= -\tfrac12\partial_x a^2$ | ponderomotive | ½ present and correct ✓ | + +--- + +## 4. Findings + +Ordered by severity. **Status: CONFIRMED** = derivation and code verified; **SUSPECTED** = probable issue, evidence stated. + +### F1 — CONFIRMED (root cause): `normalization.py` $v_0$ is $\sqrt2$ larger than the engine's velocity unit + +`normalization.py:91` ($v_0 = \sqrt{2T_0/m_e}$) vs. the engine convention $\sqrt{T_0/m_e}$ established by `helpers.py:69`, the collision operators, the dispersion tests, and all configs (§2). Since the physical temperature of the simulated plasma is $T_{phys} = m_e\,\sigma_{code}^2\,v_0^2 = T_{0,code}\cdot(m_e v_0^2)$: + +- **Under (a) as written: $T_{phys} = 2\,T_{0,code}\,T_{0,ref}$.** A `normalizing_temperature: 2000eV`, `T0: 1.0` run is a 4000 eV plasma. +- $L_0 = \sqrt2\lambda_{De}$, so all `dim="x"`/`dim="k"` string conversions carry a spurious $\sqrt2$ (box sizes, gradient scale lengths in `datamodel.py:159–177`, laser $k_0$ at `simulation.py:81`). +- The species-config `T0` parameter is effectively in units of $2T_{0,ref}$ — half the naïve expectation. + +The dimensionless dynamics of numeric-input runs are unaffected (the solver core never sees $v_0$); the damage is to the *physical interpretation* of every run and to dimensional-string inputs. + +### F2 — CONFIRMED (propagation): logged units and regression fixtures lock in the $\sqrt2$ values + +`write_units()` (`modules.py:119–141`) logs `v0`, `x0`, `box_length`, `c_light`, `nuee` computed under convention (a). For the 2000 eV resonance case the fixtures record `c_light: 11.302`, `v0: 2.652e7 m/s`, `x0: 12.14 nm` (`tests/test_vlasov1d/test_config_regression/resonance_derived_config.yml:125–132`, same in the fokker-planck and multispecies fixtures); the convention-(b) truth values are $c/v_{th} = 15.98$, $v_{th} = 1.876\times10^7$ m/s, $\lambda_{De} = 8.585$ nm. Additionally `nuee` is evaluated at $T_{0,ref}$ while the plasma actually simulated is at $2T_{0,ref}$, so the logged collisionality does not describe the simulated plasma. Any fix of F1 must regenerate these fixtures. + +### F3 — CONFIRMED: EM branch — $\hat c$ is $\sqrt2$ too small; partial cancellation hides it for wavelength drivers only + +$\beta = 1/\hat c$ with $\hat c = c/v_0$ (`modules.py:58,121`) feeds the wave solver and the CFL limit. Under (a), $\hat c = c/(\sqrt2 v_{th})$ — $\sqrt2$ smaller than the engine's thermal unit warrants. For **wavelength-specified** drivers the $\sqrt2$'s cancel in the product $\hat c\hat k = \frac{c}{\sqrt2 v_{th}}\cdot\sqrt2 k\lambda_{De}$ — which is why `test_em_dispersion.py` and `srs.yaml`'s EM branch behave physically. For **numeric** AKW drivers (`srs-debug-small.yaml`: `k0: 1.0, w0: 2.79` used as-is at `simulation.py:59`), $\hat k$ is *not* rescaled while $\hat c$ still carries the $\sqrt2$, so the EM dispersion $\hat\omega^2 = \hat n_e + \hat c^2\hat k^2$ is $\sqrt2$-inconsistent with the ES branch's $\hat k = k\lambda_{De}$ interpretation. Mixed-input SRS runs are therefore internally inconsistent between branches. + +### F4 — CONFIRMED: Krook target Maxwellian is hard-coded to $e^{-v^2/2}$ (convention-locked, $T{=}1$, electron grid only) + +`fokker_planck.py:245`. Three separate problems: (i) it bakes in convention (b) with variance exactly 1, so it must be changed **in lockstep** with any convention fix or it will spuriously heat/cool by 2× in temperature; (ii) even today, it ignores the species `T0` and `mass` — any species initialized at $T_0 \ne 1$ (e.g. `twostream.yaml`, $T_0 = 0.2$) is dragged toward $T = 1$; (iii) it is built on the electron grid and applied only to `"electron"` (`fokker_planck.py:171`). It is a fixed-target BGK operator: conserves density only, not momentum or energy. Currently `is_on: False` in shipped configs, so latent. + +### F5 — CONFIRMED: `Dougherty.compute_vbar` is missing the $1/n$ normalization + +`fokker_planck.py:81` returns $\int v f\,dv = n\bar u$ while its docstring claims "Mean velocity ⟨v⟩". `discrete_temperature` (`driftdiffusion.py:101ff`) *does* divide by $\int f\,dv$ — the two moments are asymmetric. Consequence: the drag centers on $n\bar u$, the operator relaxes toward $e^{-\beta(v - n\bar u)^2}$, and **momentum is not conserved where $n(x) \ne 1$**; $T_{target}$ is also biased by the wrong centering. Negligible for the shipped EPW/SRS configs ($n \approx 1 \pm 10^{-4}$), real for bump-on-tail / large density perturbations with collisions on. Fix: divide by $\int f\,dv$. + +### F6 — CONFIRMED: storage central moments centered on $n\bar u$, not $\bar u$; misleading names + +`storage.py:140–143`: the field moment `"v"` is $\int v f\,dv = n\bar u$ (a flux), and `p`, `q` are centered on it. Only exact when $n = 1$; biased for `nlepw-ic.yaml` (10% density perturbation) and `bump-on-tail.yaml`. Same integrand is named `"v"` in field moments but `"mean_j"` in scalars (`storage.py:304`) — one of the labels is wrong; `j` is the honest name (or divide by $n$ and keep "v"). Note this is the same class of bug as F5, appearing independently in two places. + +### F7 — CONFIRMED: the two `-flogf` entropy diagnostics have opposite signs + +Field moment `storage.py:144` computes $+\int f\log|f|\,dv$ (no minus, uses $f$); scalar `storage.py:307` computes $-\int|f|\log|f|\,dv$. For $f > 0$ these are exact negatives. The label `-flogf` matches the scalar; the field version is sign-flipped. + +### F8 — SUSPECTED (intent unclear): super-Gaussian $\alpha$ fixes $\langle v^4\rangle/\langle v^2\rangle$, not the variance, for $m \ne 2$ + +`helpers.py:72`: $\alpha = \sqrt{3\Gamma(3/m)/\Gamma(5/m)}$ normalizes the kurtosis ratio $\langle v^4\rangle/\langle v^2\rangle = 3v_{th}^2$ for all $m$, but the variance equals $v_{th}^2$ only at $m = 2$. Numerically the realized variance is $\{1.240, 1.371, 1.449\}\times T_0/m$ for $m = \{3,4,5\}$ — so the measured `p`/`n` moment will *not* equal the input `T0` for super-Gaussian species. May be an intentional flat-top-EDF temperature definition; if so it should be documented, because the `T0` config docstring ("Temperature") and the second-moment diagnostic disagree with it. All shipped configs use $m = 2$, where it is exact. + +### F9 — CONFIRMED (gap/trap): collision-frequency normalization chain is entirely manual + +- `approximate_ee_collision_frequency` returns **Hz** and is only *logged* (`modules.py:119,132`); nothing converts it to code units or feeds it to the FP operator. +- The FP/Krook rate magnitudes (`baseline`, `bump_height`) are taken as **raw floats** (`functions.py:97–98`), not passed through `normalize()` — users must supply $\nu$ already in code units ($1/\omega_{p0}$), which is nowhere documented (`config.md` lists `baseline` with no units). +- The correct conversion is $\hat\nu = \nu_{ee}[\mathrm{Hz}]\cdot\tau = \nu_{ee}/\omega_{p0}$ with **no** $2\pi$ (the NRL rate is a true s⁻¹ rate); `_tf1d/modules.py:78` has this pattern (`nuee_norm = nuee/wp0`), `_vlasov1d` has no analog. A user thinking in cyclic frequency is one step from a silent $2\pi$ error. +- The LB/Dougherty $\nu$ and the NRL $\nu_{ee}$ agree only up to an O(1) factor — identifying them silently is itself an approximation worth a docs note. +- $\nu$ is a prescribed space-time envelope; it does not track the local $n(x)/T(x)^{3/2}$ (limitation, not a bug). + +Suggested: log a `nuee_norm` alongside `nuee`, and document `baseline`'s units. + +### F10 — CONFIRMED: species `T0` and drift `v0` bypass the unit machinery + +`simulation.py:213–214` take `float(cfg.T0)`, `float(cfg.v0)`; `datamodel.py:20–21` type them as bare floats. The `normalize(dim="temp")` and `dim="v"` branches are dead code for `_vlasov1d`. Consequences: `T0: "500 eV"` is impossible (only the global `normalizing_temperature` is dimensional), and the drift `v0` is expressed in units of the reference $v_0$ while the thermal width is in $\sigma$ units — under convention (a) these differ by $\sqrt2$ *within the same distribution*, a genuine user trap (drift "1.0" is not "one thermal width"). + +### F11 — CONFIRMED (latent): `dx = xmax/nx` ignores `xmin` + +`grid.py:55` should be `(xmax − xmin)/nx`; the x-axis itself (`grid.py:74`) spans $[x_{min}, x_{max}]$ with the wrong spacing when $x_{min} \ne 0$, and the `kx` grid (`grid.py:77`) and the semi-Lagrangian interpolation period (`pushers/vlasov.py:31`, `period = xmax`) inherit the error. `modules.py:112` uses the correct $(x_{max}-x_{min})$, highlighting the discrepancy. All shipped configs use `xmin: 0.0`, so latent but real. + +### F12 — Minor items (confirmed, low impact) + +1. **`AmpereSolver` class docstring** (`field.py:230`) claims $j = \sum_s (q_s/m_s)\int vf_s\,dv$; the code (`field.py:263–264`) correctly omits $1/m_s$. Docstring bug only. +2. **Energy diagnostics lack ½ and a total**: `mean_P`, `mean_e2`, `mean_de2` are $2\times$ the respective energies (consistently, so `mean_P + mean_e2` is still conserved), and no diagnostic sums kinetic + field energy. A dedicated conservation monitor would have caught F1 earlier. +3. **`dw0 = 0.0 # ???`** placeholder at `simulation.py:87` for the intensity/wavelength driver. +4. **Dead config knob**: `GridConfig.c_light` (`datamodel.py:73`; set in `wavepacket.yaml`) is never read — `c_light`/`beta` are always recomputed from the normalization. +5. **`grid.t` axis spacing**: `nt = int(tmax/dt + 1)`, `tmax = dt·nt` overshoots the request by up to ~dt and `linspace(0, tmax, nt)` has spacing $\ne$ dt. Cosmetic (saves use their own axes). +6. **EM plasma term is electron-only** (`vector_field.py:292–301`, `field.py:152`): ions omitted from the transverse current ($\sim m_e/m_i$ — fine, but an asymmetry vs. Poisson/Ampère which include all species). +7. **Point-source amplitude** uses the vacuum $k = \omega/c$ (`field.py:67`); realized plasma amplitude larger by $k_{vac}/k_{plasma}$ — already documented in the class docstring. +8. **Half-cell axis convention**: `vax`/`x` hold cell *centers* while the extents name the *edges* ($\pm v_{max}$) — keep in mind when labeling axes. + +### Verified correct (checked because they looked suspicious, and passed) + +- The entire collisionless solver core: all coefficients unity, signs consistent (§3.3). +- The ponderomotive chain: the $-\tfrac12\partial_x a^2$, the $(q^2/m^2)$ factor, and the $c/v_0$ scaling of $a_0$ combine *exactly* — no cycle-average ½ missing, no hidden $\sqrt2$. +- The Buet "factor of 2" in the Dougherty operator: correctly absorbed by $\beta = 1/2T$; equilibrium width equals the measured variance exactly. +- The $E_x$ driver amplitude $E = \omega a_0$: consistent vector-potential convention, matches the docs. +- Midpoint-rule ($\mathrm{sum}\cdot dv$) integration: used uniformly in init, moments, and field solves. + +--- + +## 5. Recommended resolution and lockstep checklist + +**Adopt convention (b) globally**: in `electron_debye_normalization`, set $v_0 = \sqrt{T_0/m_e}$ (so $L_0 = \lambda_{De}$), and redefine `vth_norm()` accordingly. This leaves the initializer, all collision operators, all tests, and all numeric configs untouched and correct, and makes the logged physical units true. + +Things that must change together / be re-verified: + +1. `normalization.py:91–92` ($v_0$, $x_0$) — **and nothing else in that file. Do NOT redefine `vth_norm()` (:48–50)**: the codebase-wide audit (`ADEPT_CONVENTIONS_AUDIT.md`, finding C1) found that `vth_norm()` has exactly four callers, all in `vfp1d/`, which is self-consistently built on the $\sqrt{2T/m}$ convention (its initializer carries a compensating factor). Changing `vth_norm()` does nothing for `_vlasov1d` (never called here) and would silently halve VFP-1D's initialized temperature. +2. Regenerate the three `*_derived_config.yml` regression fixtures (they lock `c_light`, `v0`, `x0` at the $\sqrt2$ values). +3. Re-check every config that uses **dimensional string** inputs (`srs.yaml`: `xmax: 100um`, gradient scale length `200um`, laser wavelength/intensity) — their physical meaning shifts by $\sqrt2$ (they become *correct*; the previously-inferred physical parameters of past runs were what was wrong). +4. `c_light`/`beta` changes for all EM runs: re-validate `test_em_dispersion.py`, the point-source amplitude test, and numeric-driver SRS configs (`srs-debug-small.yaml` `w0`/`k0` values were presumably tuned under the old $\hat c$ — see F3). +5. Krook target (`fokker_planck.py:245`) — no change needed under (b), but fix its $T_0$/species handling anyway (F4) so it stops being convention-locked. +6. Docs: state the normalization explicitly in `docs/source/solvers/vlasov1d/` (velocity unit $=\sqrt{T_0/m_e}$, $L_0 = \lambda_{De}$, $\hat k = k\lambda_{De}$, Maxwellian $e^{-v^2/2}$ at $T{=}1$), and document `baseline` units (F9) and the drift-`v0` units (F10). +7. ~~Audit sibling solvers for the same $v_0$ dependency~~ **Done — see `ADEPT_CONVENTIONS_AUDIT.md`** (codebase-wide audit with per-solver immunity verdicts). Summary: `electron_debye_normalization` is consumed only by `_vlasov1d`, `_vlasov2d`, `_pic1d` — all σ-convention, all fix-safe, and only `_vlasov1d` has fixtures to regenerate. `vfp1d` genuinely runs in $\sqrt{2T/m}$ units but via `laser_normalization` + `vth_norm()`, so it is untouched by the scoped fix (see amended item 1). `_tf1d` has a private inline duplicate of the √2 (`_tf1d/modules.py:62`) that must be fixed separately. All other modules are immune. + +Independent bug fixes (any order, no convention coupling): F5 (`compute_vbar` $1/n$), F6 (moment centering + naming), F7 (`-flogf` sign), F11 (`dx` with `xmin`), F12.1 (docstring), F12.2 (add an energy-conservation scalar). + +### Suggested verification runs after any fix + +- Landau damping and ion-acoustic tests must still pass unchanged (they pin convention (b)). +- A Dougherty-collisions run with a strong density perturbation: check $\partial_t\int v f\,dv \approx 0$ (exercises F5). +- An SRS wavelength-driver run: confirm backscatter resonance moves to the physically correct location once `c_light` and the µm conversions change together. +- Compare the second-moment temperature diagnostic against `T0` for an $m = 4$ super-Gaussian to decide F8's intended convention. + +--- + +*Audit method: four parallel independent reviews (initialization/grids/derived quantities; solver core; collisions; diagnostics/docs/tests/configs), followed by cross-checking of all load-bearing claims against the source. No source files were modified.* From d3a47c9a7adb98aae397a328060dc1ee25dff2c6 Mon Sep 17 00:00:00 2001 From: Phil Travis Date: Mon, 20 Jul 2026 11:52:06 -0700 Subject: [PATCH 11/11] jk: removed audit files --- .../2026-07-16_ADEPT_CONVENTIONS_AUDIT.md | 225 ---------------- .../2026-07-16_VLASOV1D_CONVENTIONS_AUDIT.md | 249 ------------------ 2 files changed, 474 deletions(-) delete mode 100644 docs/audits/2026-07-16_ADEPT_CONVENTIONS_AUDIT.md delete mode 100644 docs/audits/2026-07-16_VLASOV1D_CONVENTIONS_AUDIT.md diff --git a/docs/audits/2026-07-16_ADEPT_CONVENTIONS_AUDIT.md b/docs/audits/2026-07-16_ADEPT_CONVENTIONS_AUDIT.md deleted file mode 100644 index 4b1ac623..00000000 --- a/docs/audits/2026-07-16_ADEPT_CONVENTIONS_AUDIT.md +++ /dev/null @@ -1,225 +0,0 @@ -# ADEPT Codebase-Wide Physics & Conventions Audit - -**Date:** 2026-07-16 -**Scope:** every solver module and shared physics file in `adept/` — `_vlasov2d`, `_pic1d`, `_tf1d`, `_lpse2d`, `_spectrax1d`, `_hermite_poisson_1d`, `vfp1d`, `osiris`, `vlasov1d2v`, plus the shared `normalization.py`, `electrostatic.py`, `functions.py`, `driftdiffusion.py`, `utils.py`, and all tests/configs. -**Companion document:** `VLASOV1D_CONVENTIONS_AUDIT.md` (the `_vlasov1d` deep audit that established the root cause). This document extends that audit to the rest of the codebase and answers: *which solvers are immune?* -**Method:** six parallel independent module audits, each cross-checked against source; load-bearing claims re-verified by hand. No source files modified. - ---- - -## 1. Executive summary - -**The √2 bug is confined to the `electron_debye_normalization` consumers — and even there, only to the unit-conversion/reporting layer.** The dimensionless *dynamics* of every solver in the codebase are internally self-consistent. The full damage inventory of `normalization.py:91` ($v_0 = \sqrt{2T_0/m_e}$, $L_0 = \sqrt2\lambda_{De}$): - -- `_vlasov1d`, `_vlasov2d`, `_pic1d`: logged physical units √2-wrong, physical temperature 2× the label, dimensional-string inputs √2-off, EM-branch $\hat c$ √2-small (details per module below). -- `_tf1d`: **independently re-implements the identical √2** in its own `write_units` (`_tf1d/modules.py:62`) — same symptom, separate code, needs its own fix. -- Everyone else: **immune** (see the table). - -**One critical amendment to the fix prescription in the vlasov1d audit:** `vth_norm()` (`normalization.py:48–50`) must **NOT** be redefined. Grep-verified, it has exactly four callers, all in `vfp1d/`, which is built self-consistently on the $v_{th} \equiv \sqrt{2T/m}$ (most-probable-speed) convention — its initializer carries a compensating factor (`vfp1d/helpers.py:111`: $\alpha = \sqrt{3\Gamma(3/m)/(2\Gamma(5/m))} = 1$ at $m{=}2$, vs. `_vlasov1d`'s $\sqrt2$) so the realized variance is exactly $T/m$. Dropping the 2 in `vth_norm()` would do nothing for `_vlasov1d` (which never calls it) and would **silently halve VFP-1D's initialized temperature** while leaving its transport coefficients at $T_0$. The safe fix touches only `electron_debye_normalization` lines 91–92. - -### 1.1 Immunity table (the answer to "which solvers are immune") - -"Bug" = the current √2 in `normalization.py:91`. "Fix" = changing `electron_debye_normalization` to $v_0=\sqrt{T_0/m_e}$ (only; `vth_norm()` untouched). - -| Solver | Working vth convention | Affected by bug NOW? | Affected by fix? | Notes | -|---|---|---|---|---| -| `_vlasov1d` | $\sqrt{T/m}$ (σ) | **YES** — logged units, physical T (2×), dimensional strings, EM $\hat c$; fixtures lock wrong values | Fix-safe for tests; must regenerate 3 fixtures | See companion doc | -| `_vlasov2d` | $\sqrt{T/m}$ (σ) | **YES (partial)** — same chain: logged units, physical T, EM $\hat c$; ES dynamics immune | **Fix-safe** — all 3 tests self-consistent, **no fixtures** | Cleanest of the affected: already fixed several vlasov1d bugs | -| `_pic1d` | $\sqrt{T/m}$ (σ) | **Logged units only** — dynamics immune (all inputs numeric) | **Fix-safe** — no fixtures, no test depends on it | Loader ≡ vlasov1d initializer convention | -| `_tf1d` | $\sqrt{T/m}$ (σ) | Not via `normalization.py` — but its **own private copy** of the √2 (`modules.py:62`) corrupts its logged units identically | Fix does **nothing** here; `modules.py:62` needs its own `2.0*T0 → T0` | Dynamics/tests fully decoupled from both | -| `_lpse2d` | $\sqrt{T/m}$ (σ), physical units | **IMMUNE** — never touches `normalization.py` | **IMMUNE** | Own ps/µm unit system; Bohm-Gross & Landau literature-exact | -| `_spectrax1d` | AW-Hermite $\alpha=\sqrt{2T/m}$ (self-consistent) | **IMMUNE** — zero imports of `normalization.py` | **IMMUNE** | α is a raw config float; physics identical to σ-convention | -| `_hermite_poisson_1d` | AW-Hermite $\alpha=\sqrt{2T/m}$ (self-consistent) | **IMMUNE** — zero imports | **IMMUNE** | Not even registered in ergoExo dispatch | -| `vfp1d` | $\sqrt{2T/m}$ (most-probable, self-consistent) | **IMMUNE** — uses `laser_normalization` ($v_0=c$), never `electron_debye_normalization` | **IMMUNE if** `vth_norm()` is left alone; **BROKEN if** `vth_norm()` is "fixed" | The reason the fix must not touch `vth_norm()` | -| `osiris` wrapper | OSIRIS native ($u_{th}=\sqrt{T/m}/c$, in-deck) | **IMMUNE** — `skin_depth_normalization` ($T_0$=None, $v_0=c$); no eV→uth conversion exists in the wrapper | **IMMUNE** | `vth_norm()` unreachable (would TypeError on $T_0$=None) | -| `vlasov1d2v` | $\sqrt{T/m}$ (σ), $m$≡1 | **IMMUNE by orphan status** — not in ergoExo dispatch, would KeyError before running | **IMMUNE** | Legacy code; carries many un-fixed bugs (§3.9) | - -**Fully immune to both bug and fix:** `_lpse2d`, `_spectrax1d`, `_hermite_poisson_1d`, `osiris`, `vfp1d` (conditional on the fix not touching `vth_norm()`), and `vlasov1d2v` (by virtue of being unrunnable). -**Affected now:** `_vlasov1d`, `_vlasov2d`, `_pic1d` (via `normalization.py`) and `_tf1d` (via its private duplicate). -**Every solver's dimensionless dynamics are unaffected in shipped configs/tests** — the corruption is confined to logged units, physical interpretation of results, dimensional-string inputs, and the EM-branch $\hat c$. - -### 1.2 Why no test ever caught the √2 - -The bug survived every test suite in the repository, and understanding why is itself a finding — the test suites, by construction, cannot see this class of error: - -1. **Every physics test works entirely in dimensionless code units.** Landau damping, ion-acoustic, Bohm-Gross, two-stream, gyro — all set $k$, $\omega$, $T_0$, box sizes as numeric floats, which pass through `normalize()` untouched (`normalization.py:57–58`). The buggy conversion layer is simply never on the tested code path. A purely numeric run *is* a correct simulation under the σ-convention; only its physical labels are wrong, and no test reads the labels. -2. **The theory references the tests compare against live in the same code-unit world.** The kinetic roots come from `electrostatic.py` with `maxwellian_convention_factor=2` and $v_{th}=1$ — i.e., the reference is expressed in the engine's internal convention, so agreement validates internal consistency, not the dimensional dictionary. The tests pin *which* convention the engine uses (invaluable for this audit) but cannot detect that `normalization.py` speaks a different one. -3. **Where a test does touch the dimensional layer, it is self-referential.** `_vlasov2d`'s `test_em_dispersion` computes $\hat c$ from the same `electron_debye_normalization` the solver uses and places the driver at $\omega^2 = 1+\hat c^2 k^2$ with that same $\hat c$ — bug and reference move in lockstep, so it passes with the wrong $c$. -4. **The only tests that pin dimensional values pin the *wrong* ones.** The `_vlasov1d` `*_derived_config.yml` regression fixtures were generated under the buggy convention, so they actively *enforce* the √2 values (`c_light: 11.302`, `x0: 12.14 nm`, …) — a fix makes tests fail, not the bug. -5. **No solver has an absolute physical-units validation** (e.g., a Landau rate checked against a rate in Hz for a stated density and temperature, or an SRS resonance checked against a wavelength in nm), and none has a total-energy conservation diagnostic. Either would have caught the mismatch the first time a dimensional input mattered. - -The general lesson: dimensionless-physics tests validate the engine; only a test that crosses the units boundary — dimensional input in, dimensional observable out, compared against an independent physical reference — can validate the normalization layer. The repo currently has zero such tests. Adding one per normalization entry point (a "round-trip units test") is the cheapest structural guard against recurrence, and belongs in the fix PR alongside item 6 of §4. - -### 1.3 The three thermal-velocity conventions in the codebase - -All are *internally* self-consistent within their modules; the mixing hazard is at module boundaries and in shared helpers: - -| Convention | Definition | Used by | -|---|---|---| -| σ (RMS / std-dev) | $v_{th} = \sqrt{T/m}$, Maxwellian $e^{-v^2/2}$ at $T{=}1$, $L_0=\lambda_{De}$ | `_vlasov1d`, `_vlasov2d`, `_pic1d`, `_tf1d`, `_lpse2d`, `vlasov1d2v` dynamics; `electrostatic.py` (mcf=2 default); OSIRIS `uth` | -| Most-probable | $v_{th} = \sqrt{2T/m}$ | `vfp1d` (with compensating init factor); `vth_norm()`; **`electron_debye_normalization` (the outlier)**; `_tf1d/modules.py:62` (private duplicate, diagnostics-only) | -| AW-Hermite scale | $\alpha = \sqrt{2T/m}$ as basis parameter; represented Maxwellian still has variance $T/m$ | `_spectrax1d`, `_hermite_poisson_1d` (Schumer–Holloway / Parker–Dellar standard) | - ---- - -## 2. Normalization wiring map - -Grep-verified callers of each `normalization.py` entry point: - -| Entry point | Definition | Callers | -|---|---|---| -| `electron_debye_normalization` | $v_0=\sqrt{2T_0/m_e}$ ← **the bug**, $L_0=v_0/\omega_{p0}$ | `_vlasov1d`, `_vlasov2d` (`modules.py:55`), `_pic1d` (`simulation.py:78`) | -| `laser_normalization` | $v_0=c$, $L_0=c/\omega_L$, $n_0=n_{crit}$, $T_0$ not self-consistent with $v_0$ (documented) | `vfp1d` (`base.py:17`) only | -| `skin_depth_normalization(_from_frequency)` | $v_0=c$, $L_0=c/\omega_{p0}$, $T_0$=None | `osiris` only | -| `vth_norm()` | $\sqrt{2T_0/m_0}/v_0$ — hard-codes most-probable | **`vfp1d` only** (grid.py:118, base.py:93, 98, 152) | -| `speed_of_light_norm()` | $c/v_0$ | `_vlasov1d`, `_vlasov2d`, `_pic1d` (√2-tainted); `osiris` (=1, exact) | -| `normalize(s, norm, dim)` | numeric passthrough; strings via $L_0$/τ/$v_0$/$T_0$ | vlasov family, `functions.py`, `vfp1d/helpers.py` profiles | -| — (no normalization import) | | `_tf1d` (own inline pint), `_lpse2d` (own unit system), `_spectrax1d`, `_hermite_poisson_1d` (raw α floats), `vlasov1d2v` (reads cfg keys nothing populates) | - -Shared physics kernels: -- `electrostatic.py` — dispersion/Z-function utilities. Self-consistent under σ-convention: `maxwellian_convention_factor=2` (default, used by **every** caller in the repo) means $f_0\propto e^{-v^2/2v_{th}^2}$, $\xi=\omega/(\sqrt2 k v_{th})$. This is the reference that pins the working convention of `_vlasov1d`, `_vlasov2d`, `_pic1d`, `_tf1d`, `_spectrax1d`, `_hermite_poisson_1d` tests. -- `driftdiffusion.py` — Dougherty/LB kernel; measures $T=\langle(v-\bar v)^2\rangle$ (or spherical $\langle v^4\rangle/3\langle v^2\rangle$), relaxes to that width. Convention-agnostic and self-adapting; verified clean in both audits. -- `functions.py` — envelope/profile layer for the vlasov family. Numeric inputs pass through; dimensional strings go through $L_0$ (√2-tainted until the fix). One bug found (§3.10). -- `utils.py` — no physics content. - ---- - -## 3. Per-module findings - -Finding IDs continue the companion doc's F-series where the bug is a copy; new IDs are per-module. - -### 3.1 `_vlasov2d` — affected now (units/EM), fix-safe, partially cleaned-up lineage - -Convention: σ, identical to `_vlasov1d` (`helpers.py:62` `v_th=√(T0/mass)`, same super-Gaussian α; 2-D init correctly normalized, $\iint f\,d^2v = n$, equal widths in $v_x,v_y$). Landau test pins σ=1 via the shared kinetic roots. Pushers all unity-coefficient; TE-mode Maxwell curl signs verified; initial 2-D Poisson correct; magnetic rotation $\theta=-(q/m)B_z dt$ verified by the gyro test. - -Bug-copy status vs `_vlasov1d`: **F5 fixed** (vbar divides by n, `fokker_planck.py:32`), **F6 fixed** (moments centered on $u=j/n$, `storage.py:100`), **F11 fixed** (`dx=(xmax-xmin)/nx`, `grid.py:60`), F7 N/A (no entropy diagnostic). **F4 still present** (Krook target $e^{-v_x^2/2}e^{-v_y^2/2}$ hard-coded, `fokker_planck.py:127–131`; latent, off in configs). **F10 still present** (`float(cfg.T0/v0x/v0y)`, `simulation.py:158–163`). **F8 present** (same α, latent at m=2). - -√2 exposure: calls `electron_debye_normalization` (`modules.py:55`) → logged units √2-wrong (F2 analog), physical T 2× label (F1 analog), EM $\hat c = c/v_0$ √2-small (F3 analog). `test_em_dispersion` passes *because it is self-consistent, not correct* — it recomputes the same wrong $\hat c$ for the driver. **No regression fixtures exist**, so the fix requires no fixture regeneration; all three tests pass unchanged (the EM test moves in lockstep with the fix). - -New (minor): `mean_KE` has the ½, field proxies $\langle E^2\rangle,\langle B^2\rangle$ don't (no combined conserved-energy diagnostic); `Txx/Tyy` are pressures ($n\cdot$variance) despite the T-name, and `T` omits the mass factor (fine for electrons); the separable per-axis Dougherty never isotropizes $T_x \leftrightarrow T_y$ (no cross-axis coupling operator exists in this module); dead config knobs `UnitsConfig.laser_wavelength`, `Z`. - -### 3.2 `_pic1d` — dynamics immune, logged units affected, fix-safe - -Convention: σ. The particle loaders (`helpers.py:34` quiet inverse-CDF, `:58` random rejection) use `v_thermal = np.sqrt(T0/mass)` — byte-for-byte the vlasov1d initializer convention. Pinned independently by `test_bohm_gross.py:34` (kinetic dielectric with vth=1) and `test_landau_damping.py:38–41` ($\omega^2=1+3k^2$, σ=1 Landau rate). - -Verified clean: deposit/gather use the identical B-spline kernel (momentum-conserving, self-force-free); uniform loading deposits exactly $n=1$ in expectation ($w=n_0L/N$, partition of unity); Poisson $E_k=-i\rho_k/k$ unity prefactor; KDK leapfrog and Yoshida4 push coefficients exact; ponderomotive $(q/m)^2(-\tfrac12\partial_x a^2)$ matches vlasov1d's verified chain. - -√2 exposure: `write_units` (`modules.py:44–57`) logs √2-wrong `v0, x0, c_light, box_length` (a 2000 eV epw run is physically 4000 eV). Dynamics immune: every shipped config and test input is numeric; loader reads raw `T0`. **No fixtures.** Fix is strictly an improvement — nothing breaks. - -New findings: -- **P1 (confirmed):** energy diagnostics mix extensivity: `mean_KE = 0.5·m·Σw v²` is a box integral with the ½; `mean_e2 = mean(e²)` is a per-cell mean without ½ or dx. `mean_KE + mean_e2` is not conserved; no valid energy monitor exists. -- **P2 (confirmed):** quiet and random loaders truncate velocity differently for drifting species — quiet clips to $[-v_{max}, v_{max}]$ absolute (`helpers.py:38`), random clips to $[v_0-v_{max}, v_0+v_{max}]$ (`helpers.py:65`). Same config, two different realized distributions; quiet wrongly clips the high-v side of a drifting Maxwellian. Latent for shipped cold-beam decks. -- **P3 (confirmed):** `mean_KE/mean_p` are sums but named `mean_*` (naming, cf. F6-class). -- Inherited: F11 (`dx=xmax/nx` while particle wrap and placement correctly use $x_{max}-x_{min}$ — the field grid and particle box disagree for $x_{min}\ne0$), F10, F8. - -### 3.3 `_tf1d` — decoupled from `normalization.py`, but carries a private copy of the √2 - -Convention: σ. Derived from the pushers: linearizing continuity + momentum ($-u\partial_x u - \frac{1}{n}\partial_x(p/m) - \frac{q}{m}E$) + adiabatic energy ($\gamma=3$) + Poisson gives $\omega^2 = 1+3k^2$, exactly what `test_resonance.py:19` asserts; the kinetic branch uses the shared roots (mcf=2). The Poisson sign ($E_k=+i\rho_k/k$, i.e. $\partial_x E=-\rho$) and the momentum sign ($-(q/m)E$) are *both* opposite the textbook and cancel — verified stable by derivation and by the passing test. Landau closure decays the field at exactly $2\,\mathrm{Im}\,\omega$, matching `test_landau_damping.py:67`. - -**T1 (confirmed, the headline):** `_tf1d/modules.py:62` re-implements `v0 = √(2·T0/m_e)` inline in its own `write_units()` (with its own pint registry — no import of `normalization.py`). Every tf1d run logs `v0, x0, c_light, beta, box_length, sim_duration` √2-wrong and implies 2× the stated temperature. Damage confined to logged diagnostics (no dimensional-string inputs exist in tf1d; the EM `WaveSolver` branch is commented out; no fixtures). **Fixing `normalization.py` does not fix this** — `modules.py:62` needs its own `2.0*T0 → T0` in lockstep. - -Other findings: -- **T2 (confirmed):** the kinetic-γ pressure closure (`pushers.py:195` `wr_corr=(wrs²-1)/k²` + γ forced to 1) yields $\omega^2 = (w_{rs}^2-1)T_0+1$ — exact only at $T_0=1$. Convention- and T₀-locked, same class as vlasov1d's Krook (F4). Latent (all configs use $T_0=1$). -- **T3 (confirmed):** `modules.py:78` `nuee_norm = nuee/wp0` — the *correct* Hz→code conversion — is computed and then **discarded** (never logged, never used). Meanwhile `physics..trapping.nuee` is an unrelated raw ML-input float. There is no collisional friction in the tf1d equations at all. (F9-class trap.) -- **T4 (confirmed):** `docs/source/usage/tf1d.md` momentum equation has a spurious $1/n$ on the E-force and mislabels the Poisson equation. Docs-only. -- **T5 (suspected, latent):** `EnergyStepper` evolves $p$ but advects/compresses $p/m$ — ambiguous mass normalization of the pressure variable; inconsistent with the momentum equation for mobile ions ($m=1836$). All configs have ions off. -- **T9 (minor):** `resonance_search.yaml` sets `ion.landau_damping: True` with `ion.is_on: False` (inert, misleading). - -### 3.4 `_lpse2d` — fully immune; thermally clean; unrelated bugs found - -Physical-units solver (ps/µm, Gaussian-cgs-derived) with its own `write_units` (`helpers.py:71`); zero imports of `normalization.py`. Convention: $v_{te} = c\sqrt{T_e/511}$ = $\sqrt{T_e/m_e}$ (σ) used consistently in the Bohm-Gross term ($e^{-i\,1.5\,v_{te}^2 k^2/\omega_{pe}\,dt}$, coefficient 3/2 correct for σ), the Landau damping rate (verified literature-exact: $\sqrt{\pi/8}(1+\tfrac32 k^2\lambda_D^2)(k\lambda_D)^{-3}e^{-3/2-1/(2k^2\lambda_D^2)}$), the sound speed, and the TPD threshold. Numerically cross-checked: the test config's driver $\omega = 1.5k^2v_{te}^2/\omega_{p0} = 19.7 \approx 20$ as configured. TPD/SRS coupling constants contain no thermal factor (convention-independent); laser $E_0=\sqrt{8\pi I/c}$ and WKB swelling $(1-n/n_c)^{-1/4}$ verified. - -Findings (none √2-class): -- **L1 (confirmed):** logged `lambda_D = vte/w0` (`helpers.py:111`) divides by the **laser** frequency instead of $\omega_{pe}$ — factor 2 too small at $n_c/4$. Diagnostic only (dynamics compute $\omega_{pe}^2/k^2v_{te}^2$ directly). -- **L2 (confirmed, latent runtime bug):** the `E2` electrostatic-driver path calls `self.epw.driver(...)` (`core/vector_field.py:84`) but the active `SpectralEPWSolver` has no `driver` attribute (only the commented-out `SpectralPotential` does) → `AttributeError` for any config with an `E2` driver. Latent only because `test_epw_frequency` is currently disabled (`pass` body) — meaning **the module's dispersion conventions have no live regression test**. -- **L3 (confirmed):** `core/trapper.py` is dead code (never instantiated by `SplitStep`); it is a stale duplicate of `_tf1d`'s trapper. Note: it uses `electrostatic.py`, *not* `driftdiffusion.py` — the vlasov1d audit's cross-module note claiming lpse2d uses the Dougherty kernel is wrong for the current tree. -- **L4 (confirmed):** the σ-convention is undocumented in the lpse2d docs — the main *future* √2 risk here is someone "harmonizing" it onto `normalization.py`. -- **L5 (minor):** the Landau/Bohm-Gross formulas are duplicated byte-for-byte in `SpectralPotential` and `SpectralEPWSolver` (drift risk); `nu_coll`'s `/2` (amplitude vs energy rate) is correct but easy to double-count. - -### 3.5 `_spectrax1d` and `_hermite_poisson_1d` — fully immune; self-consistent AW-Hermite basis - -Both use the asymmetrically-weighted Hermite basis (Schumer–Holloway / Parker–Dellar) with scale $\alpha = \sqrt{2T/m}$ taken as a **raw config float** — zero imports of `normalization.py`. The represented Maxwellian has variance $\alpha^2/2 = T/m$: physically identical to the σ-convention solvers. Verified by derivation (streaming ladder $\alpha\sqrt{n/2}$ + force ladder $\sqrt{2n}/\alpha$ + Ampère/Poisson coupling → Langmuir $\omega=\omega_{pe}$ exactly, Bohm-Gross $\omega^2=1+3(k\lambda_{De})^2$ with $\lambda_{De}=\alpha/\sqrt2$) and by tests: both modules' Landau tests set $\alpha_e = \sqrt2\,k\lambda_D/k$ explicitly against the shared mcf=2 kinetic roots (2%/5% tolerance for hermite-poisson; a cross-module test pins the two modules' E-coupling equal to 1e-12). - -The coincidence that $\alpha=\sqrt{2T/m}$ matches the `normalization.py:91` outlier is harmless — there is no code coupling in either direction. Both integrator paths (DoPri8 and Lawson-RK4 exponential operators) share identical ladder coefficients. The historical inverted E-coupling bug in hermite_poisson (`C[n+1]` vs `C[n-1]`) is fixed and regression-locked by three tests. - -Findings (all diagnostic-label/minor): **S1** logged `lambda_D` is the *total* (ion-dominated) Debye length while `k_norm` uses electron-only and hard-codes mode 1 (`base_module.py:186,189`); **S2** the ion "temperature" diagnostic is a velocity variance missing the $m_i$ factor ($T_i/m_i$, misleading next to $T_e$; `storage.py:386–405`); **S3** shipped `landau-damping.yaml` has `Nn: 4` — far too few Hermite modes to resolve the damping it is named for (tests override to 512/32); `_spectrax1d/helpers.py` is an all-stub file not on any code path. - -### 3.6 `vfp1d` — immune to the bug; the reason the fix must not touch `vth_norm()` - -Uses `laser_normalization` (`base.py:17`): $v_0 = c$, $L_0 = c/\omega_L$, $n_0 = n_{crit}$, $T_0$ = reference eV (documented as not self-consistent with $v_0$). Velocity grid in units of $c$. Thermal convention: **most-probable, $v_{th} = \sqrt{2T/m}$**, deliberately and self-consistently: - -- `vth_norm()` supplies $\sqrt{2T_0/m}/c$ (4 call sites, all vfp1d); -- the initializer's width factor `helpers.py:111` is $\alpha = \sqrt{3\Gamma(3/m)/(2\Gamma(5/m))}$ — note the extra `/2` vs the vlasov-family α — giving $\alpha=1$ at $m{=}2$ and realized variance exactly $T/m$ (the √2 and the /2 cancel); -- the v-grid extent `grid.py:118` `vmax = 8·vth_norm()/√2` = 8 standard deviations (the /√2 converts most-probable→σ); -- the IB coefficient `base.py:81` ($0.093373\,\lambda_{\mu m}^2/T_{keV}$ per $10^{15}$ W/cm²) normalizes $v_{osc}^2$ to $2T/m$, consistent; -- the temperature diagnostic `storage.py:317–322` uses the spherical variance $T = \langle v^4\rangle/3\langle v^2\rangle$ — reads the true temperature, no √2. - -Collision operators ($\nu_{ee}$ coefficient anchored to $v_0=c$ via $r_e$, Lorentz e–i with Epperlein–Haines Z*, Rosenbluth I/J integrals) are temperature-convention-independent. **No internal √2 inconsistency exists.** - -**C1 (confirmed, cross-module, high):** because items above are keyed to `vth_norm()` while the compensating α is a separate constant, redefining `vth_norm()` → $\sqrt{T_0/m}/v_0$ (as the vlasov1d audit's §5.1 originally suggested) silently initializes VFP-1D at **half** the intended temperature while all transport coefficients stay at $T_0$; the Spitzer/Epperlein–Haines gold test (`test_kappa_eh.py`, $\kappa\propto T^{5/2}$) would likely fail, and the pure-operator tests would *not* catch it. **The fix must be confined to `electron_debye_normalization:91–92`.** (The companion doc's §5.1 has been amended accordingly.) - -Other findings: **C2** `base.py:118` stores the bound method `norm.vth_norm` instead of calling it (harmless repr-string in cfg; wrong); **C3** `storage.py:170` hard-codes $9.09\times10^{21}$ cm⁻³ per $n_c$ — that is $n_{crit}(351\,\mathrm{nm})$; wrong labeling for any other `laser_wavelength` (should derive from `norm.n0`); **C4** = F11 copy (`grid.py:71` `dx=xmax/nx`); **C6/C7 (suspected, latent — IB not in any shipped config):** the production inverse-bremsstrahlung wiring looks under-normalized — `w0_norm` is identically 1.0 under laser normalization, the Langdon-factor argument $Z^2 n_i/(\omega_0 v^3)$ carries no collision-frequency coefficient, and `vosc2_per_intensity` (thermal-speed² units) is consumed as $v_{osc}^2$ in grid ($c$) units — a $\sim(c/v_{th})^2$ discrepancy between the unit tests (where the grid unit *is* the thermal speed) and production. Recommend a Spitzer-IB validation run before enabling IB in production. - -### 3.7 `osiris` wrapper — fully immune - -Uses `skin_depth_normalization(_from_frequency)` ($v_0=c$, $T_0$=None). The wrapper drives a native OSIRIS deck: **no eV→`uth` conversion exists anywhere in the wrapper** — `uth` values pass through verbatim from the deck, so the T→v √2 trap cannot occur. Verified exact: $a_0 \to E_{peak} = a_0\omega_L m_e c/e \to I = E^2\epsilon_0 c/2$ (matches $I\lambda_{\mu m}^2 = 1.37\times10^{18}a_0^2$, test-confirmed); `c_light = beta = 1.0` identically; box lengths in true skin depths; `units.yaml` deliberately omits temperature-dependent keys (consistent with $T_0$=None). Temperature diagnostics (`plots.py:1004` $T=\sum u_{th,i}^2$ in $m_ec^2$ units; `:1059` momentum variance) carry no spurious factors. `vth_norm()` is unreachable (would TypeError on $T_0$=None). - -**C8 (low, latent):** the adaptive-box feature's `reference_density = 0.25` "quarter-critical" (`density.py:75`) assumes the deck's $n_0 = n_{crit}$ (true for LPI decks with `omega0=1`); a deck normalized to a different density would silently mis-scale the density bounds (the *length* normalization stays correct). - -### 3.8 Shared `electrostatic.py` and `functions.py` - -`electrostatic.py`: self-consistent under the σ-convention throughout. $Z$ via Faddeeva, $Z' = -2(1+\xi Z)$ correct; `maxwellian_convention_factor=2` (default, used by every caller in the repo) ⇒ $f_0\propto e^{-v^2/2v_{th}^2}$, root returned as $\xi k v_{th}\sqrt{2} = \omega$. One cosmetic blemish (**T7**): the Newton `initial_root_guess` is ω-scaled ($\sqrt{\omega_p^2+3k^2v_{th}^2}$) but the root variable is ξ-scaled — converges anyway. - -`functions.py`: numeric inputs pass through; `EnvelopeFunction` amplitudes (`baseline`, `bump_height`) are raw floats (the F9 collision-frequency trap); `SineFunction.wavenumber` correctly uses `dim="k"`. **T8 (suspected):** `LinearFunction`/`ExponentialFunction` normalize `val_at_center` — a *density* — with `dim="x"` (`functions.py:159–161, 178–180`); dimensionally wrong for any string input. - -### 3.9 `vlasov1d2v` — orphaned legacy; a museum of the known bugs plus three new ones - -Not registered in the ergoExo dispatch (`_base_.py:294–333` → NotImplementedError); reads `cfg["units"]["derived"]` keys nothing populates → cannot run. Never calls `normalization.py`. Convention: σ (init `exp(-(v_x^2+v_y^2)/2T_0)`, mass≡1; the super-Gaussian machinery is commented out so `m` is silently ignored). - -Carries live copies of: **F5** (`fokker_planck.py:93,98` — vbar *and* the diffusion coefficient computed without ÷n), **F6 + F7** (`storage.py:152–158, 296`), **F10**, **F11** (`helpers.py:211`). F4 present but dead (Krook never instantiated). - -New bugs (would matter if ever revived): -- **N1 (confirmed):** the density profile is computed for every basis and folded into the ion background (`helpers.py:287`), but the line applying it to $f$ (`helpers.py:83–85`) is **commented out** — electrons always start uniform, so non-uniform configs begin with spurious net charge and no actual density perturbation in $f$. -- **N2 (confirmed):** the explicit Dougherty substeps apply ν twice (`fokker_planck.py:134–135, 141–142`: `dfdt = nu*ddx(...)` then `f + dt*nu*dfdt`) — collisional relaxation runs at $\nu^2$. -- **N4 (confirmed):** `integrator.py:234` passes the **nu_ee** envelope args when building the **nu_ei** profile (copy-paste) — the configured nu_ei profile is silently ignored. -- **N3 (design note):** without collisions the $v_y$ dimension is dynamically inert (no $v_y$ force, streaming by $v_x$ only) — the entire $n_{vy}$ grid is dead compute except through the FP operators. -- Minor: default scalars broadcast over the $v_y$ axis only (`storage.py:292–295`). - -Recommendation: either delete `vlasov1d2v` or quarantine it with a module-level comment; in its current state it is a trap for anyone who greps for reference implementations. - -### 3.10 Cross-module systemic observations - -1. **Driver amplitude conventions differ across the Vlasov family** (N5): `_vlasov1d` Ex driver $\propto \omega a_0$; `_vlasov2d` current source $\propto \omega^2 a_0$; `vlasov1d2v` $\propto |k| a_0$. Same config key, three physical meanings. Worth unifying or documenting per-solver. -2. **No solver has a valid total-energy conservation diagnostic.** vlasov1d/vlasov2d/pic1d all mix ½-factors and extensivity between kinetic and field terms. A per-solver conserved-energy scalar is the single cheapest guard that would have caught the √2 class of bugs. -3. **The F-series bugs are lineage-correlated:** `_vlasov2d` (newest refactor) fixed F5/F6/F11; `_pic1d` inherits F10/F11 by importing `_vlasov1d`'s grid/simulation; `vlasov1d2v` (oldest) has everything. When fixing `_vlasov1d`, fix the copies in the same PR: Krook targets (`_vlasov1d`, `_vlasov2d`, dead `vlasov1d2v`), `dx`/xmin (`_vlasov1d/grid.py:55`, `vfp1d/grid.py:71`, `vlasov1d2v/helpers.py:211`), raw-float T0/v0 (all three vlasov-family simulation.py files). -4. **Hz→code-unit collision-frequency conversion** exists correctly in exactly one place (`_tf1d/modules.py:78`) and is dead there; nowhere is it live. Every solver that takes a collision frequency takes it as a raw code-unit float with undocumented units. - ---- - -## 4. Consolidated fix checklist (supersedes §5 of the vlasov1d audit where they differ) - -**The √2 fix, safely scoped:** -1. `normalization.py:91–92`: $v_0 \to \sqrt{T_0/m_e}$, $x_0 \to \lambda_{De}$ — inside `electron_debye_normalization` **only**. -2. **Do NOT change `vth_norm()`** (`normalization.py:48–50`) — it is vfp1d-private and correct for vfp1d (C1). If its name is judged misleading, rename to `most_probable_speed_norm()`; do not change the value. -3. `_tf1d/modules.py:62`: `2.0*T0 → T0` (private duplicate, T1). -4. Regenerate the three `_vlasov1d` `*_derived_config.yml` regression fixtures (only vlasov1d has fixtures). -5. Fix the convention-locked Krook targets in the same PR (`_vlasov1d/…/fokker_planck.py:245`, `_vlasov2d/…/fokker_planck.py:127–131`): build from species T0/mass. -6. Re-validate: vlasov1d Landau + ion-acoustic (must pass unchanged), vlasov2d Landau/gyro/EM-dispersion (EM moves in lockstep), pic1d all three (unchanged), vfp1d `test_kappa_eh` (must pass unchanged — the sentinel that `vth_norm()` was left alone), numeric-driver SRS configs re-tuned per vlasov1d F3. -7. Document the per-solver conventions in each solver's docs page (σ vs most-probable vs AW-α; see §1.3) — the biggest residual risk is a future "harmonization" that imports the wrong convention into a currently-clean module (especially `_lpse2d`, L4). - -**Independent bug fixes, by priority:** -- High (live, wrong physics if exercised): vfp1d IB normalization audit before production use (C6/C7); lpse2d `E2` driver path + re-enable `test_epw_frequency` (L2). -- Medium (live, wrong diagnostics): lpse2d logged `lambda_D` (L1); pic1d/vlasov2d/vlasov1d energy-diagnostic normalization + add conserved-energy scalars; spectrax ion-T mass factor (S2); vfp1d hard-coded $n_{crit}(351\,nm)$ (C3); pic1d quiet-loader drift truncation (P2). -- Low / hygiene: tf1d docs (T4) and dead `nuee_norm` (T3 — log it instead of dropping it); `functions.py` `val_at_center` dim (T8); vfp1d bound-method store (C2); spectrax `Nn:4` config (S3); delete or quarantine `vlasov1d2v` and lpse2d's dead trapper (L3); unify driver-amplitude conventions or document them (N5). - ---- - -## 5. What was verified clean (highlights) - -- All solver cores' dimensionless coefficients: vlasov1d/2d pushers and field solves, pic1d push/deposit/gather/Poisson, tf1d fluid system (including its double sign cancellation, verified by derivation), Hermite ladder algebra in both spectral modules (verified by derivation to the Langmuir and Bohm-Gross limits), lpse2d Bohm-Gross/Landau/TPD/SRS coefficients against literature, vfp1d collision coefficients and temperature diagnostics, osiris a0↔intensity conversions. -- The shared `driftdiffusion.py` Dougherty kernel (both audits): self-adapting, Buet factor-2 correctly absorbed, spherical temperature correct. -- `electrostatic.py`'s Z-function algebra and mcf=2 semantics — the single reference that consistently pins six modules' tests to the σ-convention. - ---- - -*Companion deep-dive for `_vlasov1d`: `VLASOV1D_CONVENTIONS_AUDIT.md`. Method: six parallel independent module audits (one Opus agent per slice) with all load-bearing claims re-verified against source before inclusion. No source files were modified.* diff --git a/docs/audits/2026-07-16_VLASOV1D_CONVENTIONS_AUDIT.md b/docs/audits/2026-07-16_VLASOV1D_CONVENTIONS_AUDIT.md deleted file mode 100644 index 81a302a5..00000000 --- a/docs/audits/2026-07-16_VLASOV1D_CONVENTIONS_AUDIT.md +++ /dev/null @@ -1,249 +0,0 @@ -# Vlasov-1D Physics & Conventions Audit - -**Date:** 2026-07-16 -**Scope:** `adept/_vlasov1d/` (all files), plus the shared code it depends on: `adept/normalization.py`, `adept/driftdiffusion.py`, `adept/electrostatic.py`, and the `tests/test_vlasov1d/` suite and example configs. -**Trigger:** the discovery that `normalization.py:91` uses the $v_0 = \sqrt{2T_0/m}$ (most-probable-speed) convention while `_vlasov1d/helpers.py:69` uses $v_{th} = \sqrt{T_0/m}$ (RMS / standard-deviation convention). Neither has been changed yet; this audit determines which convention the module actually runs in, inventories every equation, and lists all inconsistencies found. - ---- - -## 1. Executive summary - -**The engine's working convention is $v_0 = \sqrt{T_0/m}$.** The distribution initializer, the Fokker–Planck operator, the Krook operator, the Landau-damping and ion-acoustic tests (the module's quantitative gold standards), and every numeric example config are all mutually consistent under: code velocity in units of the thermal *standard deviation* $\sqrt{T_0/m_e}$, code length unit $L_0 = \lambda_{De}$, code wavenumber $= k\lambda_{De}$, Maxwellian $\propto e^{-v^2/2}$ at $T=1$. - -**`adept/normalization.py:91–92` is the sole outlier** (`electron_debye_normalization`: $v_0 = \sqrt{2T_0/m_e}$, $L_0 = \sqrt{2}\,\lambda_{De}$). Because numeric config inputs pass through `normalize()` untouched, this does **not** corrupt the dynamics of numeric-input runs — but it means: - -- A run with `normalizing_temperature: 2000eV` and species `T0: 1.0` is physically a **4000 eV** plasma (factor 2 in temperature). -- Every dimensional **string** input converted with $L_0$ (box sizes in µm, gradient scale lengths, laser $k_0$) is off by $\sqrt{2}$. -- Every **logged** physical unit (`v0`, `x0`, `c_light`, `box_length`) is off by $\sqrt{2}$, and the regression fixtures currently lock in those wrong values. -- The EM wave speed $\hat c = c/v_0$ fed to the wave solver is $\sqrt2$ too small relative to the engine's thermal unit (with a partial cancellation for wavelength-specified drivers; see F3). - -**Recommended fix direction:** change `normalization.py` to $v_0 = \sqrt{T_0/m_e}$ (so $L_0 = \lambda_{De}$), *not* the initializer. Fixing `helpers.py:69` instead (to $\sqrt{T_0/2m}$) would break the currently-passing Landau-damping and ion-acoustic tests and invalidate every existing config's driver $(k_0, \omega_0)$ values. Section 5 lists everything that must move in lockstep. - -The collisionless solver core (`vector_field.py`, pushers) was verified **clean**: every coefficient is exactly unity under the declared normalization and is convention-independent. The convention bug lives entirely in the dimensional-conversion layer. - -Beyond the $\sqrt2$ issue, the audit found several independent bugs and traps, listed in Section 4 (notably: FP `compute_vbar` missing $1/n$; storage central moments centered on $n u$ instead of $u$; a sign inconsistency between the two `-flogf` entropy diagnostics; the super-Gaussian $\alpha$ fixing the wrong moment for $m\ne2$; the Krook target hard-coded to $T=1$; and the collision-frequency normalization chain being entirely manual). - ---- - -## 2. The two candidate conventions - -With $\omega_{p0} = \sqrt{n_0 e^2/(\epsilon_0 m_e)}$, $\tau = 1/\omega_{p0}$ in both cases: - -| Quantity | (a) `normalization.py` as written | (b) engine's actual convention | -|---|---|---| -| Velocity unit $v_0$ | $\sqrt{2T_0/m_e}$ (most-probable speed) | $\sqrt{T_0/m_e}$ (RMS / std-dev) | -| Length unit $L_0 = v_0/\omega_{p0}$ | $\sqrt{2}\,\lambda_{De}$ | $\lambda_{De}$ | -| Maxwellian at $\hat T = 1$ | $\propto e^{-v^2}$ (variance ½) | $\propto e^{-v^2/2}$ (variance 1) | -| Code wavenumber $\hat k$ | $\sqrt2\, k\lambda_{De}$ | $k\lambda_{De}$ | -| Bohm–Gross | $\hat\omega^2 = 1 + \tfrac{3}{2}\hat k^2$ | $\hat\omega^2 = 1 + 3\hat k^2$ | -| $\hat c = c/v_0$ (2000 eV) | 11.30 | 15.98 | - -**Why the test suite never caught this:** every physics test (Landau, ion-acoustic, EM dispersion) works in dimensionless code units — numeric inputs bypass `normalize()` entirely (`normalization.py:57–58`), so the buggy conversion layer is never on the tested path, and the theory references (`electrostatic.py`, mcf=2, $v_{th}=1$) are expressed in the engine's internal convention, so agreement validates internal consistency rather than the physical-units dictionary. Worse, the only tests that *do* pin dimensional values — the `*_derived_config.yml` regression fixtures — were generated under the buggy convention and therefore lock the wrong values in. No test crosses the units boundary (dimensional input in, dimensional observable out, against an independent physical reference), which is precisely the kind of test that must accompany the fix (see `ADEPT_CONVENTIONS_AUDIT.md` §1.2 for the full analysis). - -Evidence pinning convention (b) as the working one, strongest first: - -1. **Landau-damping test** (`tests/test_vlasov1d/test_landau_damping.py:24`) calls `electrostatic.get_roots_to_electrostatic_dispersion(1.0, 1.0, k)` — i.e. $\omega_{pe}=1$, $v_{th}=1$ with `maxwellian_convention_factor=2` (`adept/electrostatic.py:79,115,98`), which is the kinetic dielectric for $f_0 \propto e^{-v^2/2v_{th}^2}$, $\xi = \omega/(\sqrt2 k v_{th})$. The measured $\omega(k)$ matches these roots to 2 decimals only if the simulated Maxwellian has $\sigma = 1$. (E.g. $k=0.3 \Rightarrow \omega = 1.1598$, which is exactly the `w0` in `resonance.yaml` and `epw.yaml`; under convention (a) the resonance would sit near 1.07 and the test would fail.) -2. **Ion-acoustic test** (`test_ion_acoustic_wave.py`) uses $c_s^2 = ZT_e/m_i$ and $\omega^2 = k^2c_s^2/(1+k^2\lambda_D^2)$ with $\lambda_D = 1$ in code units — asserting $L_0 = \lambda_{De}$. -3. **Initializer** (`helpers.py:69–76`), **Krook target** (`fokker_planck.py:245`), and **Dougherty stationary state** (`fokker_planck.py:51,96` + `driftdiffusion.py`) all build/assume $e^{-v^2/(2T/m)}$, variance $T/m$. -4. Every numeric example config (`epw.yaml`, `resonance.yaml`, `wavepacket.yaml`, `bump-on-tail.yaml`) uses driver $(k_0,\omega_0)$ pairs consistent only with $\hat k = k\lambda_{De}$ and $\sigma = \sqrt{T/m}$. - ---- - -## 3. Equation inventory - -### 3.1 Normalization layer (`adept/normalization.py`) - -| Location | Expression | Meaning | Convention note | -|---|---|---|---| -| :88 | $\omega_{p0} = \sqrt{n_0e^2/\epsilon_0 m_e}$, $\tau = 1/\omega_{p0}$ | time unit | standard | -| **:91** | $v_0 = \sqrt{2T_0/m_e}$ | velocity unit | **outlier — convention (a)** | -| :92 | $x_0 = v_0/\omega_{p0} = \sqrt2\lambda_{De}$ | length unit | inherits $\sqrt2$ | -| :48–50 | `vth_norm()` $= \sqrt{2T_0/m_0}/v_0 = 1$ | "thermal velocity" | convention (a); **never called** in `_vlasov1d` | -| :52–53 | `speed_of_light_norm()` $= c/v_0$ | $\hat c$ | inherits $\sqrt2$ (F3) | -| :37–38 | NRL $\log\Lambda_{ee}$ | Coulomb log | uses reference $T_0$ | -| :45 | $\nu_{ee} = 2.91\times10^{-6}\, n_{cc}\log\Lambda\, T_{eV}^{-3/2}$ Hz | NRL e–e rate | **diagnostic only**, never converted to code units (F9) | -| :56–74 | `normalize()`: x/L0, t/τ, v/v0, T/T0, k·L0 | dim→code | numeric inputs pass through untouched; `dim="temp"` and `dim="v"` branches are dead code for `_vlasov1d` | - -### 3.2 Initialization & grids (`helpers.py`, `grid.py`, `modules.py`, `simulation.py`, `datamodel.py`) - -| Location | Expression | Meaning | Convention note | -|---|---|---|---| -| helpers.py:65–66 | `dv = 2 vmax/nv`; cell-centered `vax` on $[-v_{max}+dv/2,\ v_{max}-dv/2]$ | v-grid | edges at $\pm v_{max}$; half-cell center offset | -| **helpers.py:69** | $v_{th} = \sqrt{T_0/m}$ | thermal width | **convention (b)** | -| helpers.py:72 | $\alpha = \sqrt{3\,\Gamma(3/m)/\Gamma(5/m)}$ ($=\sqrt2$ at $m{=}2$) | super-Gaussian width factor | fixes $\langle v^4\rangle/\langle v^2\rangle = 3v_{th}^2$; variance $=v_{th}^2$ **only** at $m=2$ (F8) | -| helpers.py:74–76 | $f \propto \exp[-|(v-v_d)/(\alpha v_{th})|^m]$ | init EDF; $m{=}2$: $e^{-(v-v_d)^2/(2T_0/m)}$ | variance $T_0/m$ — convention (b) | -| helpers.py:80,84 | `f /= sum(f)·dv`; `f *= n_prof` | $\int f\,dv = n(x)$ | midpoint rule, consistent with all moments | -| grid.py:55 | `dx = xmax/nx` | cell width | **ignores `xmin`** (F11) | -| grid.py:59–60 | `dt = min(dt, 0.95·dx/c)`, $c = 1/\beta$ | EM CFL | $\hat c$ inherits F1's $\sqrt2$ | -| grid.py:65–66,75 | `nt = int(tmax/dt+1)`; `t = linspace(0, dt·nt, nt)` | time axis | overshoot ≤ dt; `grid.t` spacing ≠ dt (F12) | -| grid.py:74,77 | cell-centered x; `kx = 2π·fftfreq(nx, d=dx)` | spectral grid | inherits F11's dx for `xmin≠0` | -| modules.py:112 | `box_length = (xmax−xmin)·L0` | box size in µm | uses `xmin` correctly (unlike grid.py:55); value carries F1's $\sqrt2$ | -| simulation.py:213–214 | `v0 = float(cfg.v0)`; `T0 = float(cfg.T0)` | species drift & temperature | bare code floats, bypass `normalize()` (F10) | -| simulation.py:74–76 | $a_0 = a_{0,std}\cdot(c/v_0)$ | intensity → quiver velocity in $v_0$ units | internally consistent with pusher (§3.3) | -| simulation.py:79–85 | $\hat k_0 = k_{phys}L_0$, $\hat\omega_0 = \omega_{phys}\tau$ | wavelength driver → code units | $\hat k_0$ carries F1's $\sqrt2$; $\hat\omega_0/\hat k_0 = \hat c$ ✓ | -| simulation.py:87 | `dw0 = 0.0 # ???` | frequency offset placeholder | unresolved (F12) | - -### 3.3 Solver core (`solvers/vector_field.py`, `solvers/pushers/vlasov.py`, `solvers/pushers/field.py`) — **verified clean** - -The dimensionless system implemented, with every coefficient exactly unity under the declared normalization (and independent of the (a)/(b) choice, since $v_0$ cancels): - -$$\partial_t f_s + v\,\partial_x f_s + \frac{\hat q_s}{\hat m_s}\left(E - \frac{\hat q_s}{2 \hat m_s}\partial_x a^2\right)\partial_v f_s = C[f_s]$$ -$$\partial_x E = \sum_s \hat q_s \int f_s\,dv \ (+\ \text{static ion background}), \qquad \partial_t E = -\sum_s \hat q_s \int v f_s\,dv$$ -$$\partial_t^2 a = \hat c^2\,\partial_x^2 a - n_e\, a + S$$ - -| Location | Code | Check | -|---|---|---| -| vlasov.py:210–214, 220 | exact spectral shift $f(x - v\,dt)$ | coefficient $v_0\tau/L_0 = 1$ ✓ | -| vlasov.py:83–87, 127–129 | `force = q·e + (q²/m)·pond; accel = force/m` | $\hat q/\hat m$ E-push; ponderomotive $(\hat q^2/\hat m^2)$ ✓ | -| vector_field.py:411 | `pond = −0.5·∂ₓ(a²)` | exact instantaneous $-\frac{1}{2m}\partial_x p_\perp^2$ with $p_\perp = qA$; the $c/v_0$ scaling of $a_0$ (simulation.py:76) makes $a$ the quiver velocity in $v_0$ units — **no hidden $\sqrt2$; verified exact** | -| field.py:203–224 | $E_k = -i\rho_k/k$ | Poisson prefactor $e^2n_0/(\epsilon_0 m_e\omega_{p0}^2) = 1$ ✓ | -| field.py:263–264, 280 | $E^{n+1} = E^n - dt\,j$ | Ampère prefactor 1 ✓ | -| field.py:337–343 | $\Delta E_k = -i\frac{q}{k}\int dv\, f_k(e^{-ikv\,dt}-1)$ | exact (Hamiltonian) Ampère along free streaming ✓ | -| field.py:146–153 | leapfrog wave eq., plasma term $-n_e a$ | dispersion $\hat\omega^2 = \hat n_e + \hat c^2\hat k^2$ ✓; $n_e$ is electron-only (~$m_e/m_i$ approx., F12) | -| vector_field.py:292–301, 341 | $n_e = -\,\hat q_e \int f_e\,dv$, time-centered | sign keyed to electron charge $-1$ ✓ | -| field.py:21–26 | $E_x$ driver $= (\omega_0{+}\delta\omega)\,a_0\sin(k_0x - \omega t)$ | vector-potential amplitude convention ($E = -\partial_t A$); matches docs (config.md:312) ✓ | -| field.py:67–69, 78–80 | point source $F_0 = 2\omega\hat c\,a_0$ | reproduces amplitude $a_0$ in vacuum; plasma value larger by $k_{vac}/k_{plasma}$, documented in docstring ✓ | - -Charge-sign conventions were checked across Poisson / Ampère / wave-equation and the force term: consistent, no compensating double-error. Current density correctly has **no** mass factor in code (but see docstring bug F12). - -### 3.4 Collisions (`solvers/pushers/fokker_planck.py`, `adept/driftdiffusion.py`) - -Operator: $\partial_t f = \nu\,\partial_v[(v - \bar v)f + T\,\partial_v f]$ (Lenard–Bernstein/Dougherty, Buet notation $\beta = 1/2T$, $D = 1/2\beta = T$, drift $C = 2\beta D(v-\bar v) = (v-\bar v)$). - -| Location | Code | Check | -|---|---|---| -| driftdiffusion.py:101ff (`discrete_temperature`) | $T = \sum f(v-\bar v)^2 dv \,/\, \sum f\,dv$ | full 2nd central moment, no factor 2, no $m$ — variance convention, matches init ✓ | -| driftdiffusion.py:170,180 | $\beta_{init} = 1/2T$; $f_{mx} = e^{-\beta(v-\bar v)^2}$ | stationary state $e^{-(v-\bar v)^2/2T}$ ✓ | -| driftdiffusion.py:333–334 | $C = 2\beta D\,(v_{edge} - \bar v)$ | $2\beta D = 1$ exactly ✓ | -| driftdiffusion.py (flux/Chang–Cooper, implicit solve) | conservative central & positivity-preserving fluxes; $(I - dt\,\nu L)f^{n+1} = f^n$ | numerics only; $\nu\,dt$ dimensionless ✓ | -| driftdiffusion.py:20–25 | Buet "extra factor of 2" note | correctly absorbed into $\beta = 1/2T$; **no stray factor survives** (verified) | -| **fokker_planck.py:81** | `compute_vbar` $= \sum f\,v\,dv$ | returns $n\bar u$, **not** $\bar u$ — missing $1/n$ (F5) | -| **fokker_planck.py:245** | Krook target $f_{mx} \propto e^{-v^2/2}$ | hard-coded variance 1 = ($T{=}1$, $m{=}1$, convention (b)) — convention-locked (F4) | -| fokker_planck.py:262–266 | $f \to f e^{-\nu dt} + n(x) f_{mx}(1 - e^{-\nu dt})$ | BGK; conserves $n$, not momentum/energy (by design) | - -Key property: the **Dougherty operator is self-adapting** — it measures $T = \langle(v-\bar v)^2\rangle$ from $f$ and relaxes toward exactly that width, so it preserves whatever convention the initializer used and needs no change under a convention fix. The **Krook operator does not** — its width is a compile-time constant (F4). - -### 3.5 Diagnostics & storage (`storage.py`) - -| Location | Code | Meaning | Note | -|---|---|---|---| -| :134–135 | $\int(\cdot)\,dv$ = `sum·dv` | midpoint rule | consistent with init normalization ✓ | -| :139 | `n` $= \int f\,dv$ | density | ✓ | -| :140 | `v` $= \int v f\,dv$ | **raw first moment $= n\bar u$**, labeled "v" | flux, not velocity (F6) | -| :141–142 | `p` $= \int (v - n\bar u)^2 f\,dv$ | "pressure" | centered on $n\bar u$, not $\bar u$ (F6); for $n{=}1$ Maxwellian at $T_0$: `p` $= T_0$ under convention (b) ✓, would read $2T_0$ under (a) | -| :143 | `q` $= \int (v-n\bar u)^3 f\,dv$ | heat flux | same centering issue | -| :144 | `-flogf` $= \int f\log|f|\,dv$ | "entropy" | **missing minus sign** (F7) | -| :307 | `mean_-flogf` $= \langle\int -|f|\log|f|\,dv\rangle$ | entropy | has the minus; opposite sign to :144 (F7) | -| :303–306, 308 | `mean_P` $=\langle\int v^2f\rangle$, `mean_j` $=\langle\int vf\rangle$, `mean_n`, `mean_q` $=\langle\int v^3 f\rangle$, `mean_f2` | raw moments | `mean_j` and field "v" are the same integrand under different names (F6) | -| :311–312 | `mean_de2` $=\langle de^2\rangle$, `mean_e2` $=\langle e^2\rangle$ | field energy proxies | no ½; kinetic `mean_P` also lacks ½, consistently — internal conservation OK, but no combined energy monitor exists (F12) | -| :154, :313 | `pond` $= -\tfrac12\partial_x a^2$ | ponderomotive | ½ present and correct ✓ | - ---- - -## 4. Findings - -Ordered by severity. **Status: CONFIRMED** = derivation and code verified; **SUSPECTED** = probable issue, evidence stated. - -### F1 — CONFIRMED (root cause): `normalization.py` $v_0$ is $\sqrt2$ larger than the engine's velocity unit - -`normalization.py:91` ($v_0 = \sqrt{2T_0/m_e}$) vs. the engine convention $\sqrt{T_0/m_e}$ established by `helpers.py:69`, the collision operators, the dispersion tests, and all configs (§2). Since the physical temperature of the simulated plasma is $T_{phys} = m_e\,\sigma_{code}^2\,v_0^2 = T_{0,code}\cdot(m_e v_0^2)$: - -- **Under (a) as written: $T_{phys} = 2\,T_{0,code}\,T_{0,ref}$.** A `normalizing_temperature: 2000eV`, `T0: 1.0` run is a 4000 eV plasma. -- $L_0 = \sqrt2\lambda_{De}$, so all `dim="x"`/`dim="k"` string conversions carry a spurious $\sqrt2$ (box sizes, gradient scale lengths in `datamodel.py:159–177`, laser $k_0$ at `simulation.py:81`). -- The species-config `T0` parameter is effectively in units of $2T_{0,ref}$ — half the naïve expectation. - -The dimensionless dynamics of numeric-input runs are unaffected (the solver core never sees $v_0$); the damage is to the *physical interpretation* of every run and to dimensional-string inputs. - -### F2 — CONFIRMED (propagation): logged units and regression fixtures lock in the $\sqrt2$ values - -`write_units()` (`modules.py:119–141`) logs `v0`, `x0`, `box_length`, `c_light`, `nuee` computed under convention (a). For the 2000 eV resonance case the fixtures record `c_light: 11.302`, `v0: 2.652e7 m/s`, `x0: 12.14 nm` (`tests/test_vlasov1d/test_config_regression/resonance_derived_config.yml:125–132`, same in the fokker-planck and multispecies fixtures); the convention-(b) truth values are $c/v_{th} = 15.98$, $v_{th} = 1.876\times10^7$ m/s, $\lambda_{De} = 8.585$ nm. Additionally `nuee` is evaluated at $T_{0,ref}$ while the plasma actually simulated is at $2T_{0,ref}$, so the logged collisionality does not describe the simulated plasma. Any fix of F1 must regenerate these fixtures. - -### F3 — CONFIRMED: EM branch — $\hat c$ is $\sqrt2$ too small; partial cancellation hides it for wavelength drivers only - -$\beta = 1/\hat c$ with $\hat c = c/v_0$ (`modules.py:58,121`) feeds the wave solver and the CFL limit. Under (a), $\hat c = c/(\sqrt2 v_{th})$ — $\sqrt2$ smaller than the engine's thermal unit warrants. For **wavelength-specified** drivers the $\sqrt2$'s cancel in the product $\hat c\hat k = \frac{c}{\sqrt2 v_{th}}\cdot\sqrt2 k\lambda_{De}$ — which is why `test_em_dispersion.py` and `srs.yaml`'s EM branch behave physically. For **numeric** AKW drivers (`srs-debug-small.yaml`: `k0: 1.0, w0: 2.79` used as-is at `simulation.py:59`), $\hat k$ is *not* rescaled while $\hat c$ still carries the $\sqrt2$, so the EM dispersion $\hat\omega^2 = \hat n_e + \hat c^2\hat k^2$ is $\sqrt2$-inconsistent with the ES branch's $\hat k = k\lambda_{De}$ interpretation. Mixed-input SRS runs are therefore internally inconsistent between branches. - -### F4 — CONFIRMED: Krook target Maxwellian is hard-coded to $e^{-v^2/2}$ (convention-locked, $T{=}1$, electron grid only) - -`fokker_planck.py:245`. Three separate problems: (i) it bakes in convention (b) with variance exactly 1, so it must be changed **in lockstep** with any convention fix or it will spuriously heat/cool by 2× in temperature; (ii) even today, it ignores the species `T0` and `mass` — any species initialized at $T_0 \ne 1$ (e.g. `twostream.yaml`, $T_0 = 0.2$) is dragged toward $T = 1$; (iii) it is built on the electron grid and applied only to `"electron"` (`fokker_planck.py:171`). It is a fixed-target BGK operator: conserves density only, not momentum or energy. Currently `is_on: False` in shipped configs, so latent. - -### F5 — CONFIRMED: `Dougherty.compute_vbar` is missing the $1/n$ normalization - -`fokker_planck.py:81` returns $\int v f\,dv = n\bar u$ while its docstring claims "Mean velocity ⟨v⟩". `discrete_temperature` (`driftdiffusion.py:101ff`) *does* divide by $\int f\,dv$ — the two moments are asymmetric. Consequence: the drag centers on $n\bar u$, the operator relaxes toward $e^{-\beta(v - n\bar u)^2}$, and **momentum is not conserved where $n(x) \ne 1$**; $T_{target}$ is also biased by the wrong centering. Negligible for the shipped EPW/SRS configs ($n \approx 1 \pm 10^{-4}$), real for bump-on-tail / large density perturbations with collisions on. Fix: divide by $\int f\,dv$. - -### F6 — CONFIRMED: storage central moments centered on $n\bar u$, not $\bar u$; misleading names - -`storage.py:140–143`: the field moment `"v"` is $\int v f\,dv = n\bar u$ (a flux), and `p`, `q` are centered on it. Only exact when $n = 1$; biased for `nlepw-ic.yaml` (10% density perturbation) and `bump-on-tail.yaml`. Same integrand is named `"v"` in field moments but `"mean_j"` in scalars (`storage.py:304`) — one of the labels is wrong; `j` is the honest name (or divide by $n$ and keep "v"). Note this is the same class of bug as F5, appearing independently in two places. - -### F7 — CONFIRMED: the two `-flogf` entropy diagnostics have opposite signs - -Field moment `storage.py:144` computes $+\int f\log|f|\,dv$ (no minus, uses $f$); scalar `storage.py:307` computes $-\int|f|\log|f|\,dv$. For $f > 0$ these are exact negatives. The label `-flogf` matches the scalar; the field version is sign-flipped. - -### F8 — SUSPECTED (intent unclear): super-Gaussian $\alpha$ fixes $\langle v^4\rangle/\langle v^2\rangle$, not the variance, for $m \ne 2$ - -`helpers.py:72`: $\alpha = \sqrt{3\Gamma(3/m)/\Gamma(5/m)}$ normalizes the kurtosis ratio $\langle v^4\rangle/\langle v^2\rangle = 3v_{th}^2$ for all $m$, but the variance equals $v_{th}^2$ only at $m = 2$. Numerically the realized variance is $\{1.240, 1.371, 1.449\}\times T_0/m$ for $m = \{3,4,5\}$ — so the measured `p`/`n` moment will *not* equal the input `T0` for super-Gaussian species. May be an intentional flat-top-EDF temperature definition; if so it should be documented, because the `T0` config docstring ("Temperature") and the second-moment diagnostic disagree with it. All shipped configs use $m = 2$, where it is exact. - -### F9 — CONFIRMED (gap/trap): collision-frequency normalization chain is entirely manual - -- `approximate_ee_collision_frequency` returns **Hz** and is only *logged* (`modules.py:119,132`); nothing converts it to code units or feeds it to the FP operator. -- The FP/Krook rate magnitudes (`baseline`, `bump_height`) are taken as **raw floats** (`functions.py:97–98`), not passed through `normalize()` — users must supply $\nu$ already in code units ($1/\omega_{p0}$), which is nowhere documented (`config.md` lists `baseline` with no units). -- The correct conversion is $\hat\nu = \nu_{ee}[\mathrm{Hz}]\cdot\tau = \nu_{ee}/\omega_{p0}$ with **no** $2\pi$ (the NRL rate is a true s⁻¹ rate); `_tf1d/modules.py:78` has this pattern (`nuee_norm = nuee/wp0`), `_vlasov1d` has no analog. A user thinking in cyclic frequency is one step from a silent $2\pi$ error. -- The LB/Dougherty $\nu$ and the NRL $\nu_{ee}$ agree only up to an O(1) factor — identifying them silently is itself an approximation worth a docs note. -- $\nu$ is a prescribed space-time envelope; it does not track the local $n(x)/T(x)^{3/2}$ (limitation, not a bug). - -Suggested: log a `nuee_norm` alongside `nuee`, and document `baseline`'s units. - -### F10 — CONFIRMED: species `T0` and drift `v0` bypass the unit machinery - -`simulation.py:213–214` take `float(cfg.T0)`, `float(cfg.v0)`; `datamodel.py:20–21` type them as bare floats. The `normalize(dim="temp")` and `dim="v"` branches are dead code for `_vlasov1d`. Consequences: `T0: "500 eV"` is impossible (only the global `normalizing_temperature` is dimensional), and the drift `v0` is expressed in units of the reference $v_0$ while the thermal width is in $\sigma$ units — under convention (a) these differ by $\sqrt2$ *within the same distribution*, a genuine user trap (drift "1.0" is not "one thermal width"). - -### F11 — CONFIRMED (latent): `dx = xmax/nx` ignores `xmin` - -`grid.py:55` should be `(xmax − xmin)/nx`; the x-axis itself (`grid.py:74`) spans $[x_{min}, x_{max}]$ with the wrong spacing when $x_{min} \ne 0$, and the `kx` grid (`grid.py:77`) and the semi-Lagrangian interpolation period (`pushers/vlasov.py:31`, `period = xmax`) inherit the error. `modules.py:112` uses the correct $(x_{max}-x_{min})$, highlighting the discrepancy. All shipped configs use `xmin: 0.0`, so latent but real. - -### F12 — Minor items (confirmed, low impact) - -1. **`AmpereSolver` class docstring** (`field.py:230`) claims $j = \sum_s (q_s/m_s)\int vf_s\,dv$; the code (`field.py:263–264`) correctly omits $1/m_s$. Docstring bug only. -2. **Energy diagnostics lack ½ and a total**: `mean_P`, `mean_e2`, `mean_de2` are $2\times$ the respective energies (consistently, so `mean_P + mean_e2` is still conserved), and no diagnostic sums kinetic + field energy. A dedicated conservation monitor would have caught F1 earlier. -3. **`dw0 = 0.0 # ???`** placeholder at `simulation.py:87` for the intensity/wavelength driver. -4. **Dead config knob**: `GridConfig.c_light` (`datamodel.py:73`; set in `wavepacket.yaml`) is never read — `c_light`/`beta` are always recomputed from the normalization. -5. **`grid.t` axis spacing**: `nt = int(tmax/dt + 1)`, `tmax = dt·nt` overshoots the request by up to ~dt and `linspace(0, tmax, nt)` has spacing $\ne$ dt. Cosmetic (saves use their own axes). -6. **EM plasma term is electron-only** (`vector_field.py:292–301`, `field.py:152`): ions omitted from the transverse current ($\sim m_e/m_i$ — fine, but an asymmetry vs. Poisson/Ampère which include all species). -7. **Point-source amplitude** uses the vacuum $k = \omega/c$ (`field.py:67`); realized plasma amplitude larger by $k_{vac}/k_{plasma}$ — already documented in the class docstring. -8. **Half-cell axis convention**: `vax`/`x` hold cell *centers* while the extents name the *edges* ($\pm v_{max}$) — keep in mind when labeling axes. - -### Verified correct (checked because they looked suspicious, and passed) - -- The entire collisionless solver core: all coefficients unity, signs consistent (§3.3). -- The ponderomotive chain: the $-\tfrac12\partial_x a^2$, the $(q^2/m^2)$ factor, and the $c/v_0$ scaling of $a_0$ combine *exactly* — no cycle-average ½ missing, no hidden $\sqrt2$. -- The Buet "factor of 2" in the Dougherty operator: correctly absorbed by $\beta = 1/2T$; equilibrium width equals the measured variance exactly. -- The $E_x$ driver amplitude $E = \omega a_0$: consistent vector-potential convention, matches the docs. -- Midpoint-rule ($\mathrm{sum}\cdot dv$) integration: used uniformly in init, moments, and field solves. - ---- - -## 5. Recommended resolution and lockstep checklist - -**Adopt convention (b) globally**: in `electron_debye_normalization`, set $v_0 = \sqrt{T_0/m_e}$ (so $L_0 = \lambda_{De}$), and redefine `vth_norm()` accordingly. This leaves the initializer, all collision operators, all tests, and all numeric configs untouched and correct, and makes the logged physical units true. - -Things that must change together / be re-verified: - -1. `normalization.py:91–92` ($v_0$, $x_0$) — **and nothing else in that file. Do NOT redefine `vth_norm()` (:48–50)**: the codebase-wide audit (`ADEPT_CONVENTIONS_AUDIT.md`, finding C1) found that `vth_norm()` has exactly four callers, all in `vfp1d/`, which is self-consistently built on the $\sqrt{2T/m}$ convention (its initializer carries a compensating factor). Changing `vth_norm()` does nothing for `_vlasov1d` (never called here) and would silently halve VFP-1D's initialized temperature. -2. Regenerate the three `*_derived_config.yml` regression fixtures (they lock `c_light`, `v0`, `x0` at the $\sqrt2$ values). -3. Re-check every config that uses **dimensional string** inputs (`srs.yaml`: `xmax: 100um`, gradient scale length `200um`, laser wavelength/intensity) — their physical meaning shifts by $\sqrt2$ (they become *correct*; the previously-inferred physical parameters of past runs were what was wrong). -4. `c_light`/`beta` changes for all EM runs: re-validate `test_em_dispersion.py`, the point-source amplitude test, and numeric-driver SRS configs (`srs-debug-small.yaml` `w0`/`k0` values were presumably tuned under the old $\hat c$ — see F3). -5. Krook target (`fokker_planck.py:245`) — no change needed under (b), but fix its $T_0$/species handling anyway (F4) so it stops being convention-locked. -6. Docs: state the normalization explicitly in `docs/source/solvers/vlasov1d/` (velocity unit $=\sqrt{T_0/m_e}$, $L_0 = \lambda_{De}$, $\hat k = k\lambda_{De}$, Maxwellian $e^{-v^2/2}$ at $T{=}1$), and document `baseline` units (F9) and the drift-`v0` units (F10). -7. ~~Audit sibling solvers for the same $v_0$ dependency~~ **Done — see `ADEPT_CONVENTIONS_AUDIT.md`** (codebase-wide audit with per-solver immunity verdicts). Summary: `electron_debye_normalization` is consumed only by `_vlasov1d`, `_vlasov2d`, `_pic1d` — all σ-convention, all fix-safe, and only `_vlasov1d` has fixtures to regenerate. `vfp1d` genuinely runs in $\sqrt{2T/m}$ units but via `laser_normalization` + `vth_norm()`, so it is untouched by the scoped fix (see amended item 1). `_tf1d` has a private inline duplicate of the √2 (`_tf1d/modules.py:62`) that must be fixed separately. All other modules are immune. - -Independent bug fixes (any order, no convention coupling): F5 (`compute_vbar` $1/n$), F6 (moment centering + naming), F7 (`-flogf` sign), F11 (`dx` with `xmin`), F12.1 (docstring), F12.2 (add an energy-conservation scalar). - -### Suggested verification runs after any fix - -- Landau damping and ion-acoustic tests must still pass unchanged (they pin convention (b)). -- A Dougherty-collisions run with a strong density perturbation: check $\partial_t\int v f\,dv \approx 0$ (exercises F5). -- An SRS wavelength-driver run: confirm backscatter resonance moves to the physically correct location once `c_light` and the µm conversions change together. -- Compare the second-moment temperature diagnostic against `T0` for an $m = 4$ super-Gaussian to decide F8's intended convention. - ---- - -*Audit method: four parallel independent reviews (initialization/grids/derived quantities; solver core; collisions; diagnostics/docs/tests/configs), followed by cross-checking of all load-bearing claims against the source. No source files were modified.*