fix: Vlasov1d vth conventions

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

adept/_vlasov1d/datamodel.py

@@ -71,7 +71,6 @@ class GridConfig(BaseModel):
     vmax: float | None = None
     vmin: float | None = None
     parallel: tuple[str, ...] | bool = False
-    c_light: float | None = None
 
     @field_validator("parallel", mode="before")
     @classmethod

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:

adept/_vlasov1d/helpers.py

@@ -72,7 +72,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
+    # <v^4>/<v^2> = 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))

adept/_vlasov1d/modules.py

@@ -118,6 +118,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()
 
@@ -131,6 +134,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,
@@ -238,10 +242,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

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.
     """
 

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 <v>, shape (...)
+            vbar: Mean velocity <v> = ∫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
 

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, :]

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

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,17 +78,21 @@ 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)
 
     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)

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:

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

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
@@ -535,7 +557,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 |

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

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

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