Repository navigation
MaxwellWeakAmpere: StiffnessPreconditioner as default; document the discretization of MaxwellWeakAmpere and VlasovAmpereCoupling - #819
Merged
Conversation
…ling Semi-discrete systems in struphy operator notation (M1, M2, C, Lambda^1, W), Crank-Nicolson update, SchurSolver blocks A, B, C, resulting Schur system, accumulated matrix/vector and conserved quantities. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…w docs] The Schur operator M1 + dt^2/4 C^T M2 C is preconditioned by StiffnessPreconditioner(curl, mid=dt^2/4 M2, stab=M1), built per time step in prepare_time_step. Iterations hardly depend on dt (hollow cylinder 16x16x4, p=3: 39 instead of 411 PCG iterations at dt = 100 h); for dt <~ h the mass preconditioner stays slightly better. algo="explicit" falls back to the mass preconditioner, solver="cg" drops the preconditioner with a warning. New test file test_maxwell_weak_ampere.py. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Contributor
|
📖 Docs preview: https://struphy-hub.github.io/docs-preview/pr-819/ |
…nt preconditioners and solvers [preview docs] - remove the iteration-count sentence from the MaxwellWeakAmpere docstring - full-path cross-references (Derham.curl, SchurSolver, prepare_time_step, kernel functions in their *_kernels modules); feectools inverse as literal - dev-feec.rst: new sections Preconditioners (struphy.feec.preconditioner, struphy.polar.kron_2dx1d, multigrid preconditioner) and Linear solvers (schur_solver, solver) - dev-pic.rst: vlasov_maxwell accumulation kernel, new Pusher kernels section with push_v_with_efield Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The raw docstrings of SchurSolver, SchurSolverFull and SchurSolverFull3 used \\text and \\right (a LaTeX line break in a raw string), which broke the math; rewritten with pmatrix. Markdown links replaced by plain names, and the SchurSolver Parameters section now matches the signature (precond, solver_params). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…n [preview docs] Lynch, Rice and Thomas (1964) for the method, Sangalli and Tani (2016) for preconditioning spline stiffness matrices; short description of the generalized eigendecomposition used by KroneckerSumSolver. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…nalization [preview docs] Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Contributor
|
📖 Docs preview: https://struphy-hub.github.io/docs-preview/pr-819/ |
spossann
added this pull request to stack #822
October 23, 2026 08:58
spossann
added a commit
that referenced
this pull request
Oct 23, 2026
… with quasi-neutral limit; model VlasovMaxwellQuasiNeutral (#821) 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 #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 #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](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
MaxwellWeakAmpereandVlasovAmpereCouplingnow describe the discretization in struphy operator notation: semi-discrete systems (SchurSolverblocksvlasov_maxwellkernel and the conserved energies.MaxwellWeakAmpere:precond="StiffnessPreconditioner"(new literalLiteralOptions.OptsMaxwellPrecond). The Schur operatorStiffnessPreconditioner(mass_ops, "curl", sigma=1, mid=dt²/4 M2, stab="M1")(PolarStiffnessPreconditionerwith polar splines), built for the current dt inprepare_time_stepand rebuilt when dt changes (same pattern asShearAlfvenPropagator).algo="explicit"falls back to the mass preconditioner ofsolver="cg"(which takes no preconditioner, and crashed before with any preconditioner) now drops the preconditioner with a warning.src/struphy/propagators/tests/test_maxwell_weak_ampere.py: same solution and exact energy conservation with stiffness and mass preconditioner, fewer iterations with stiffness at dt ≫ h, rebuild on dt change, polar type, explicit fallback, options.PCG iterations for the Schur system (16×16×4 cells, p=3, tol 1e-8, h = 1/16, c = 1):
For dt ≲ h (e.g. CFL-limited PIC runs such as the VlasovMaxwellOneSpecies examples, dt/h ≈ 0.3) the two are about equal; the mass preconditioner remains available and is recommended there in the docstring.
VlasovAmpereCouplingkeeps the mass preconditioner: its Schur operator has no derivatives (a mass matrix with weight $1 + \Delta t^2\alpha^2 n/(4\varepsilon^2)$).Docs
MaxwellWeakAmpere/VlasovAmpereCouplingdocstrings resolve (full paths, e.g.Derham.curl,SchurSolver, the kernels in their*_kernelsmodules); the iteration-count sentence was removed from theMaxwellWeakAmperedocstring.dev-feec.rst: new sections Preconditioners (struphy.feec.preconditioner,struphy.polar.kron_2dx1d, multigrid preconditioner) and Linear solvers (schur_solver,solver). This also fixes the previously broken preconditioner links innumerics-geomFE.rstand inShearAlfvenPropagator,CurlCurlSolve,ImplicitDiffusion.dev-pic.rst: thevlasov_maxwellaccumulation kernel and a new Pusher kernels section (push_v_with_efield).SchurSolver,SchurSolverFull,SchurSolverFull3: the math did not render (\\text/\\rightin raw strings is a LaTeX line break); rewritten withpmatrix, and theSchurSolverparameters now match the signature.StiffnessPreconditioner: references for fast diagonalization (Lynch, Rice & Thomas 1964; Sangalli & Tani 2016) and the formula of the inverse.doc/gets parsed by Sphinx (crashKeyError: 'anchorname'); keep it outsidedoc/.Known issue (pre-existing, not addressed here)
The implicit
MaxwellWeakAmperestep with polar splines fails inSchurSolver.__call__(_schur += BC, domain mismatch inSumLinearOperator), also with the oldPolarMassMatrixPreconditioner. The polar test only checks that the polar stiffness preconditioner is built.Tests run locally
test_maxwell_weak_ampere.py: serial (6 passed) andmpirun -n 2(5 passed, 1 mpi_skip)test_verif_Maxwell.py,test_verif_VlasovMaxwellOneSpecies.py(serial)-m single): Maxwell, ColdPlasma (2 ranks); ColdPlasmaVlasov, LinearVlasovMaxwellOneSpecies, VlasovMaxwellOneSpecies (4 ranks)-n): no unresolved references from the changed docstringsDocs preview triggered with
[preview docs].🤖 Generated with Claude Code