Skip to content

Particles5D: draw mu with the local |B_0|, B0="equil" default (fixes #841) - #878

Merged
spossann merged 5 commits into
develfrom
particles5d-local-B0
Nov 3, 2026
Merged

spossann merged 5 commits into
develfrom
particles5d-local-B0

Conversation

@spossann

@spossann spossann commented Nov 3, 2026 •

Copy link
Copy Markdown
Member

Fixes #841.

Particles5D drew the magnetic moment as mu = -log(1-r) * vth_perp**2 / B0 with the constant LoadingParameters.B0 (default 2.0). In a tokamak the markers were therefore not Maxwellian in v_perp at their position, which added δf noise. Meanwhile the sampling density s0 already accepted a callable B0, so passing one made the draw and the weights inconsistent.

What changes

  • B0 option: LoadingParameters.B0 and GyroMaxwellian2D(B0=...) take a positive number (also NumPy/CuPy scalars), a callable of stacked logical positions (3, N), or "equil". "equil" is the new default for both: it is the local |B_0| of the species' equilibrium.
  • Loading (draw_markers): mu is drawn with B0 at the final marker positions. The disc inversion and the pole clamp now run before velocity sampling. A specific marker with a prescribed position but a drawn mu keeps its drawn v_perp.
  • One B0 for the draw and s0: Particles5D.loading_B0 is sampling_density.eval_B0, so f_init/s0 doesn't depend on mu when the moments match.
  • KineticBackground.attach_equil (forwarded by sums and scalar multiples): Particles5D attaches its equilibrium to the background, the initial condition and the sampling density. Attaching a different equilibrium to a background that already has one (one object shared by two simulations) logs a warning. The class-level default background is deep-copied per species.
  • Evaluation: GyroMaxwellian2D.eval_B0 caches the last positions, so f0 and velocity_jacobian_det at the same markers cost one |B_0| evaluation. On a phase-space meshgrid, B0 is evaluated on the spatial points only and broadcast. A constant B0 is not evaluated at all.
  • Examples: the GuidingCenter tokamak and screw-pinch examples and the LinearMHDDriftkineticCC lowmode_tae example dropped their constant-B0 workarounds (B0=B_min, B0=1.0, B0=equil.absB0) and use the default. The stale # B0=equil.absB0 comments in the cyclone, ITG and diocotron examples are gone.
  • Usage guide: a pitfall entry in the setup-simulation skill.

Why the numbers change (regression baselines)

The baselines of GuidingCenter, DriftKineticElectrostaticAdiabatic, DriftKineticMaxwell, LinearMHDDriftkineticCC and ToyDrift were updated for 1 and 2 ranks. The default params use HomogenSlab, whose |B_0| = 1 replaces the old constant B0 = 2.0. So mu (and sum_w_v2) roughly doubles, and the hybrid's mu-dependent current coupling changes accordingly. All other model baselines are unchanged (all 20 model regressions pass on 1 and 2 ranks).

Behaviour changes to note

  • A GyroMaxwellian2D with the default B0="equil" that is evaluated on its own (e.g. background.plot(...) before the simulation, outside a Particles5D species) raises a clear ValueError. Pass a number or callable there.
  • Cost: a δf run with a GyroMaxwellian2D background now evaluates the analytic |B_0| once per weight update. Measured with 2e5 markers, serial, ad hoc torus, markers moving: update_weights took 32 ms with a constant B0 and 174 ms with "equil"; it was 349 ms before the cache. Runs that already passed B0=equil.absB0 paid this before. A possible follow-up is to evaluate the projected absB0_h with the marker kernels instead.
  • mu is drawn with the analytic |B_0|, the same as f_init, while the pushers use the projected absB0_h. The difference is the usual O(h^p) projection error.

Tests

  • New src/struphy/pic/tests/test_loading_B0.py, on 1 and 4 ranks:
    • mu |B_0(x)| / vth**2 has mean 1 (Exp(1)) for markers both below and above the median field strength in a hollow torus, for pseudo-random/uniform and Sobol/disc loading. With a constant B0 the means are off.
    • With matching moments the weights are |det DF| / s_eta / Np up to 1e-10; a mismatch gives velocity-dependent weights.
    • Validation of the B0 options, scalar broadcast and meshgrid broadcast.
  • Existing tests: kinetic_background, test_output, show_distribution, scalars, splitting order, binning (serial); GC pushers, 5D current coupling, DriftKineticMaxwell coupling, pole (4 ranks); all 20 model regressions (1 and 2 ranks); the 14 affected examples with --test-run (2 ranks).
  • test_show_distribution: the (v_par, mu) bins now reach mu = 5 (mean mu = 1 with |B_0| = 1), and the noise-dominated (eta1, eta2) tolerance went from 0.31 to 0.33. That error was already 0.308 with B0 = 2.0.

🤖 Generated with Claude Code

LoadingParameters.B0 and GyroMaxwellian2D(B0=...) accept a number, a
callable of stacked logical positions (3, N) or "equil" (new default for
both), the local field strength |B_0| of the species' equilibrium.

- draw_markers draws mu = v_perp^2/(2 B0) with B0 at the final marker
  positions (disc inversion and pole clamp now before velocity sampling;
  specific markers with a prescribed position keep their drawn v_perp)
- Particles5D.loading_B0 and the sampling density s0 share one B0
  (GyroMaxwellian2D.eval_B0), so draw and weights are consistent
- KineticBackground.attach_equil (forwarded by sums/scalar multiples);
  Particles5D attaches its equilibrium to background and initial condition
- regression baselines of the 5D models updated (default HomogenSlab has
  |B_0| = 1 instead of the old constant 2.0)
- examples: GC tokamak/screw pinch and DKCC lowmode_tae use the default
  instead of constant B0 workarounds

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🟡 Changes recommended

Equilibrium attachment mutates shared distributions, causing inconsistent particle weights across populations, and the B0 cache returns stale field values after callable replacement.

2 open issues
What changed in this PR

This PR addresses #841 by drawing Particles5D magnetic moments using the local equilibrium field strength, keeping marker loading consistent with the sampling density.

Changes:

  • Adds constant, callable, and default "equil" field options, with equilibrium attachment and cached evaluation.
  • Updates marker-loading order, adds validation and statistical tests, and refreshes five model regression baselines.
  • Removes example workarounds and documents the new defaults.
File Description
src/​struphy/​pic/​tests/​test_show_distribution.py Adjusts magnetic-moment bins and statistical tolerance.
src/​struphy/​pic/​tests/​test_loading_B0.py Tests local-field loading, weights, validation, and broadcasting.
src/​struphy/​pic/​particles.py Connects distributions and loading to the species equilibrium.
src/​struphy/​pic/​base.py Finalizes positions before velocity sampling and rescales prescribed markers.
src/​struphy/​particles/​parameters.py Defines and validates the expanded B0 options.
src/​struphy/​models/​tests/​regression/​baselines/​ToyDrift.json Updates magnetic-moment regression values.
src/​struphy/​models/​tests/​regression/​baselines/​LinearMHDDriftkineticCC.json Updates energies, coupled fields, and marker baselines.
src/​struphy/​models/​tests/​regression/​baselines/​GuidingCenter.json Updates perpendicular and total energy baselines.
src/​struphy/​models/​tests/​regression/​baselines/​DriftKineticMaxwell.json Updates energy, field, and marker baselines.
src/​struphy/​models/​tests/​regression/​baselines/​DriftKineticElectrostaticAdiabatic.json Updates magnetic-moment and energy baselines.
src/​struphy/​kinetic_background/​tests/​test_maxwellians.py Supplies explicit B0 for a standalone distribution test.
src/​struphy/​kinetic_background/​maxwellians.py Adds equilibrium binding, caching, and spatial meshgrid evaluation.
src/​struphy/​kinetic_background/​base.py Adds shared field helpers and forwards equilibrium attachment through compositions.
src/​struphy/​agent_guides/​skills/​setup-simulation/​SKILL.md Documents local-field defaults and matching loading options.
examples/​ToyDrift/​diocotron_instability/​params_diocotron.py Removes stale B0 comments only.
examples/​LinearMHDDriftkineticCC/​lowmode_tae/​params_lowmode_tae.py Replaces constant-field loading with the equilibrium default.
examples/​GuidingCenter/​tokamak_equilibrium/​params_tokamak_equilibrium.py Removes explicit loading and background field overrides.
examples/​GuidingCenter/​screw_pinch_equilibrium/​params_screw_pinch_equilibrium.py Uses equilibrium defaults for loading and background.
examples/​DriftKineticElectrostaticAdiabatic/​itg_cylindre/​params_drift_kinetic.py Removes stale B0 comments only.
examples/​DriftKineticElectrostaticAdiabatic/​cyclone/​params_cyclone.py Removes stale B0 comments only.

Review effort: Balanced

Comment thread src/struphy/pic/particles.py Resolved
spossann and others added 4 commits November 3, 2026 08:59
…he key

Review of #878:
- KineticBackground.with_equil(equil, memo) replaces attach_equil: it returns
  a bound copy (GyroMaxwellian2D with B0="equil", and sums/scalar multiples
  of such) instead of changing the given object, so a background shared by
  species with different equilibria stays consistent per species. One
  memoizing with_equil in the base class (memo keeps the originals alive),
  subclasses implement _bind_equil; aliases are kept (initial condition =
  background, shared summands).
- Particles binds background and initial condition in a _bind_backgrounds
  hook right after setting them in __init__ (Particles5D: to B0_equil), so
  nothing sees unbound objects; the sampling density is bound too.
- PICVariable refers to the bound copies of its particles after allocation,
  so var.backgrounds can be evaluated after the run; the given objects stay
  unbound.
- GyroMaxwellian2D.eval_B0: the cache key includes the B0 option and the
  bound equilibrium (by identity), so replacing params["B0"] invalidates it.
- Tests: shared background with two equilibria, bound sums, cache
  invalidation, var.backgrounds after allocation.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…ached

Reverts the copy-on-bind design of 15612b6 (with_equil, _bind_equil,
_bind_backgrounds hook, PICVariable re-pointing): within a simulation there
is only one equilibrium, and attaching keeps the user's background the object
the particles evaluate. Kept from the review:
- GyroMaxwellian2D.eval_B0 cache key includes the B0 option and the attached
  equilibrium, so replacing params["B0"] invalidates it;
- attach_equil logs a warning when a different equilibrium is attached to the
  same background (shared by two simulations set up before running), instead
  of silently switching.
Tests: cache invalidation, warning, var.backgrounds attached after allocation.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@spossann
spossann merged commit 68fd601 into devel Nov 3, 2026
30 checks passed
@spossann
spossann deleted the particles5d-local-B0 branch November 3, 2026 09:31
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

In Particles5D, LoadingParameters.B0 cannot be a callable

2 participants