[spossann]

Follow-up of #819 (merged); based on devel.

Normalization (latest change)

VlasovMaxwellQuasiNeutral uses the cyclotron velocity scale ($\hat v = \hat\Omega_c \hat x$, $\hat t = 1/\hat\Omega_c$, $\varepsilon = 1$) with two new/redefined species equation parameters, passable at model initialization:

$$\lambda = \frac{\hat\Omega_c \hat x}{c},,\qquad \kappa = \frac{\hat\Omega_p \hat x}{c} = \alpha\lambda,,\qquad \lambda^2 \partial_t \mathbf E = \nabla\times\mathbf B - \kappa^2 \mathbf J,,$$

with no coefficients in the Vlasov equation and Faraday's law; the quasi-neutral limit is lam=0. (kappa was $\hat\Omega_p \hat t$ before, unused; in the light scale both definitions agree.) For a general velocity scale the coefficients are $\varepsilon^2\lambda^2$, $\kappa^2\varepsilon$, $1/\varepsilon$ and Gauss $\kappa^2/(\varepsilon\lambda^2) = \alpha^2/\varepsilon$; they reduce to $1, \alpha^2/\varepsilon, 1/\varepsilon$ for the light scale. Documented under Equation parameters in models-normalization.rst (which also had a spurious $2\pi$ in the cyclotron velocity, removed). doc_normalization is now directly below doc_pde (rule added to models/AGENTS.md).

Consequences in this PR: the propagator option kappa (displacement current) is renamed lam (coefficient $\lambda^2$), defaults $\lambda = \varepsilon\lambda_s$, $c_2 = \kappa_s^2\varepsilon$, $c_3 = 1/\varepsilon$; the model's scalars en_E, en_B, en_f are the parts of the conserved energy and en_tot is the conserved energy (en_H removed); the initial Poisson solve uses Gauss' law $\lambda^2 \nabla\cdot\mathbf E = c_2 \rho$ of the coupling step and is skipped for $\lambda = 0$ (gauss_error NaN). Example and verification test use lam, kappa=1, epsilon=-1, which are the old equations: the verification growth rate is unchanged (0.02770). Signed cyclotron units for a negative bulk charge: issue #825.

The sections below describe the propagator in its own coefficients; where they say kappa for the displacement-current coefficient, read lam (with $\lambda^2$), and en_H is now en_tot.

Summary

New propagator VlasovMaxwellCoupling: the steps of MaxwellWeakAmpere and VlasovAmpereCoupling in one implicit Crank–Nicolson step, for

$$\kappa, \mathbb M_1 \dot{\mathbf e} = c_1, \mathbb C^\top \mathbb M_2 \mathbf b - c_2, \mathbf j(\mathbf V),\qquad \dot{\mathbf b} = -\mathbb C \mathbf e,\qquad \dot{\mathbf V} = c_3, \overline{DF^{-\top}} \mathbf\Lambda^1 \mathbf e$$

(marker positions fixed), with options kappa (default 1), c1 (default 1), c2, c3 (default $\alpha^2/\varepsilon$, $1/\varepsilon$ of the species). $\kappa = 0$ drops the displacement current (quasi-neutral limit, $c_1 \nabla\times\mathbf B = c_2 \mathbf J$).

• Schur system for the mid-point field: $\big(\kappa \mathbb M_1 + \tfrac{\Delta t^2}{4}(c_1 \mathbb C^\top \mathbb M_2 \mathbb C + c_2 c_3 \mathbb A)\big), \mathbf e^{\mathrm{mid}} = \kappa \mathbb M_1 \mathbf e^n + \tfrac{\Delta t}{2}(c_1 \mathbb C^\top \mathbb M_2 \mathbf b^n - c_2 \mathbf j^n)$ with the accumulated $\mathbb A$, $\mathbf j$ of the vlasov_maxwell kernel. Then $\mathbf b^{n+1} = \mathbf b^n - \Delta t, \mathbb C \mathbf e^{\mathrm{mid}}$, velocities pushed with $\mathbf e^{\mathrm{mid}}$, $\mathbf e^{n+1} = 2\mathbf e^{\mathrm{mid}} - \mathbf e^n$ for $\kappa > 0$ (exactly CN) and $\mathbf e^{n+1} = \mathbf e^{\mathrm{mid}}$ for $\kappa = 0$ (avoids the undamped even–odd oscillation of the implicit mid-point rule for the algebraic variable $\mathbf e$).
• Well-posed for $\kappa \to 0$: the curl-curl term stays; its kernel (discrete gradients, harmonic fields) is absorbed by the plasma term $c_2 c_3 \mathbb A \approx c_2 c_3 \mathbb M_{1,n}$. Needs positive weights: with the control variate, $\kappa = 0$ is rejected (assert), $\kappa > 0$ warns. At $\kappa = 0$ the mid-point current is weakly divergence-free ($\mathbb G^\top \mathbf j^{\mathrm{mid}} = 0$); $\mathbf E$ is its Lagrange multiplier.
• Exact energy conservation of $H = \frac\kappa2 \mathbf e^\top \mathbb M_1 \mathbf e + \frac{c_1}{2} \mathbf b^\top \mathbb M_2 \mathbf b + \frac{c_2}{2c_3} \sum_p w_p |\mathbf v_p|^2$ (method energy()), for full-f weights.
• Preconditioner (default): StiffnessPreconditioner(mass_ops, "curl", sigma=kappa + dt^2 c2 c3 n_ref/4, mid=dt^2 c1/4 M2, stab="M1"), built in prepare_time_step and rebuilt when dt changes (as in MaxwellWeakAmpere: StiffnessPreconditioner as default; document the discretization of MaxwellWeakAmpere and VlasovAmpereCoupling #819). The stabilization stays positive at $\kappa = 0$. n_ref defaults to the density average $\mathbf 1^\top \mathbb A \mathbf 1 / \mathbf 1^\top \mathbb M_1 \mathbf 1$.
• Matrix-free option matrix_free (default False), as in the other 6D couplings of Matrix-free particle matrices (MarkerOperator) in the 6D coupling propagators #817: the particle matrix in the Schur operator is a MarkerOperator (vlasov_maxwell_mf + marker_operator_symm_dot), with the fillings $w_p G^{-1}$ stored per marker.
• Initial projection at κ = 0 (option initial_projection, default True): the constraint holds only at the mid-point, so its residual $\mathbf r^n = c_1 \mathbb C^\top \mathbb M_2 \mathbf b^n - c_2 \mathbf j(\mathbf V^n)$ obeys $\mathbf r^{n+1} = -\mathbf r^n$; initial data that violate $c_1 \nabla\times\mathbf B = c_2 \mathbf J$ (e.g. marker noise in $\mathbf J$) would oscillate with period 2 without damping. In the setup (first prepare_time_step, before the initial data are saved; not for a restored restart state), $(\mathbf b, \mathbf V)$ is replaced by the mid-point state of one step (implicit Euler of dt/2), which satisfies the constraint exactly and does not depend on dt; $\mathbf e$ (initial Poisson) is kept. The particle pushers of the splitting still change $\mathbf r$ by $O(\Delta t)$ per step.
• No polar splines yet (asserted).

New model VlasovMaxwellQuasiNeutral: subclass of VlasovMaxwellOneSpecies with the splitting coupling → push_eta → push_vxb and the extra scalar en_H. VlasovMaxwellOneSpecies got a _link_propagator_variables hook and a default-params hook that keys on the last propagator and the class name (output for VlasovMaxwellOneSpecies unchanged).

Verification test models/tests/verification/test_verif_VlasovMaxwellQuasiNeutral.py: Weibel instability at κ = 0 (1e5 full-f markers with sobol_antithetic loading, dt = 1, t = 160); the growth rate fitted to en_B for 80 < t < 150 must be within 7 % of the analytical 0.02784. Fitted: 0.02770 / 0.02784 / 0.02782 on 1 / 2 / 4 ranks (~45 s serial, ~20–27 s on 2 or 4 ranks). With pseudo-random markers the fit is noise-limited (about 10 % low, ~7 % seed-to-seed scatter), hence the antithetic Sobol loading.

Example examples/VlasovMaxwellQuasiNeutral/weibel_quasi_neutral (--kappa 1|0, dt = 1, i.e. 20× the dt of the VlasovMaxwellOneSpecies Weibel example).

Results

PCG iterations of one step (8×4×1 cells, p = (2,2,1), 1280 full-f markers, tol 1e-13; energy error ≤ 3e-16 relative in every case):

dt	kappa	Stiffness	Mass
0.1	1	13	10
0.1	0	20	72
2	1	19	58
2	0	22	77
20	1	23	80
20	0	23	80

$\kappa \to 0$: $|\mathbf b - \mathbf b_0| + |\mathbf V - \mathbf V_0|$ against the $\kappa = 0$ solution decreases linearly (8.4, 1.9, 0.22, 0.023 for $\kappa = 10^{-1}, \dots, 10^{-4}$).

Weibel example with sobol_antithetic loading (current version; 10⁵ full-f markers, 32 cells, t = 400, serial, ~2 min per run): γ = 0.02805 for κ = 1 and 0.02758 for κ = 0 (analytical 0.02784), max rel. change of en_H 2e-9 / 6e-9.

Earlier runs with pseudo-random markers:

run	fitted γ	max rel. change of en_H
κ = 1, dt = 1	0.0249	1.4e-9
κ = 0, dt = 1	0.0249	9.0e-10
κ = 0, dt = 0.25	0.0253	2.4e-10
κ = 0, dt = 1, with initial projection	0.0251	3.1e-9

The first three κ = 0 rows are without the initial projection (it was added later): there en_B shows a period-2 oscillation band (≈ 2e-7) for t < 80, which the projection removes (smooth en_B from t = 0, en_E 10–20× lower).

Analytical γ = 0.02784 for both κ (κγ² ≈ 8e-4 ≪ k² + 1). With pseudo-random markers the fitted rates were 10 % low, hardly dependent on dt: marker noise, removed by the antithetic Sobol loading above.

Known limitations / follow-ups

• Existing code, found in review: with clones, sobol_standard/sobol_antithetic loadings have no clone offset in the Sobol index, so all clones draw the same markers (effective sample Np / nclones). Not used by the tests here; follow-up issue.

• Strict quasi-neutrality is not enforced over the full split step: the coupling step makes $\mathbf J(\mathbf V^{\mathrm{mid}})$ divergence-free at frozen positions, but PushEta then moves the markers with $\mathbf V^{n+1}$.

• Polar splines not supported (tensor-product V1 buffers; VlasovAmpereCoupling has the same limitation).

• With the control variate (κ > 0 only) the energy is not exactly conserved and n_ref falls back to 1.

Tests run locally

(Updated after the matrix-free option and the initial projection. The branch has devel merged via #819, which also got devel merged, so the diff here is only this PR.)

• test_vlasov_maxwell_coupling.py + test_matrix_free_6d.py: serial 40 passed, mpirun -n 4 39 passed + 1 mpi_skip. New: matrix-free vs assembled steps (Colella, κ = 1 and 0, incl. n_ref), discrete equations with both variants, test_initial_projection (sign flip of r without; r = 0, dt independence, e unchanged and stationary steps with the projection). test_maxwell_weak_ampere.py on the merged state: serial and 2 ranks.
• Separate review agents for the matrix-free and the projection commits: no bugs in the matrix-free path; the projection review found that writing e_mid (an impulse of size 1/dt) into e replaced the Poisson field in the saved t = 0 data, fixed (only b and V are projected).

Earlier (first version):

• src/struphy/propagators/tests/test_vlasov_maxwell_coupling.py: serial (15 passed) and mpirun -n 4 (14 passed, 1 mpi_skip): exact energy conservation (κ = 1, 0; both preconditioners; species and explicit c1 ≠ |c2| ≠ |c3|), residual of the discrete CN Ampère and Faraday laws (with the pushed V), weakly divergence-free mid-point current at κ = 0, O(dt²) agreement with MaxwellWeakAmpere + VlasovAmpereCoupling, linear convergence in κ with bounded iterations, rebuild on dt change, option validation, control-variate rejection. A mutation check (pusher constant c2 instead of c3) makes 6 of these fail.
• Default-params model tests -m single for VlasovMaxwellQuasiNeutral and VlasovMaxwellOneSpecies (serial and 4 ranks); test_examples.py for the new example; doc/check_models_doc.py, doc/check_propagator_doc.py; ruff.
• Separate review agent on the branch: no correctness bugs; its low-severity findings (polar splines, test weaknesses) are fixed in the second commit.

🤖 Generated with Claude Code