VlasovMaxwellCoupling: implicit energy-conserving Vlasov-Maxwell step with quasi-neutral limit; model VlasovMaxwellQuasiNeutral

doc/sections/subsections/models-all.rst

@@ -21,6 +21,11 @@ Full-orbit Vlasov
     :undoc-members:
     :show-inheritance:
 
+.. autoclass:: struphy.models.VlasovMaxwellQuasiNeutral
+    :members: doc_pde, doc_scalar_quantities, doc_normalization, doc_discretization, doc_long_description, doc_examples, doc_use_cases, doc_cannot_be_used_for
+    :undoc-members:
+    :show-inheritance:
+
 .. autoclass:: struphy.models.LinearVlasovMaxwellOneSpecies
     :members: doc_pde, doc_scalar_quantities, doc_normalization, doc_discretization, doc_long_description, doc_examples, doc_use_cases, doc_cannot_be_used_for
     :undoc-members:

doc/sections/subsections/models-normalization.rst

@@ -49,7 +49,9 @@ There are four possibilities:
 
 .. math::
 
-    \hat v = v_\textnormal{c, bulk} := \hat x \Omega_\textnormal{c, bulk}/(2\pi) = \hat x\, q_\textnormal{bulk} \hat B /(m_\textnormal{bulk}2\pi)\,.
+    \hat v = v_\textnormal{c, bulk} := \hat x\, \Omega_\textnormal{c, bulk} = \hat x\, q_\textnormal{bulk} \hat B /m_\textnormal{bulk}\,,
+
+so that :math:`\hat t = 1/\Omega_\textnormal{c, bulk}`.
 
 4. ``thermal`` velocity of the bulk species,
 
@@ -106,7 +108,8 @@ In Struphy models, the following equation parameters appear:
 
 .. math::
 
-        \alpha_\textrm{s} = \frac{\hat \Omega_\textnormal{ps}}{\hat \Omega_\textnormal{cs}}\,,\qquad \varepsilon_\textrm{s} = \frac{1}{\hat \Omega_\textnormal{cs} \hat t} \,,
+        \alpha_\textrm{s} = \frac{\hat \Omega_\textnormal{ps}}{\hat \Omega_\textnormal{cs}}\,,\qquad \varepsilon_\textrm{s} = \frac{1}{\hat \Omega_\textnormal{cs} \hat t} \,,\qquad
+        \lambda_\textrm{s} = \frac{\hat \Omega_\textnormal{cs} \hat x}{c}\,,\qquad \kappa_\textrm{s} = \frac{\hat \Omega_\textnormal{ps} \hat x}{c} = \alpha_\textrm{s} \lambda_\textrm{s}\,,
 
 featuring the plasma- and cyclotron frequency of species :math:`\textrm{s}`, respectively,
 
@@ -115,6 +118,7 @@ featuring the plasma- and cyclotron frequency of species :math:`\textrm{s}`, res
         \hat\Omega_\textnormal{ps} = \sqrt{\frac{\hat n (Z_\textrm{s}e)^2}{\epsilon_0 (A_\textrm{s} m_\textnormal{H})}} \,,\qquad \hat \Omega_{\textnormal{cs}} = \frac{(Z_\textrm{s}e) \hat B}{(A_\textrm{s} m_\textnormal{H})}\,,
 
 where :math:`Z_\textrm{s}` and :math:`A_\textrm{s}` stand for the species' charge and mass number, respectively.
+
 These equation parameters are defined in :class:`~struphy.models.species.Species.EquationParameters` and can be overridden
 in the model constructor.
 

doc/sections/subsections/propagators-avail.rst

@@ -197,6 +197,11 @@ Vlasov-Maxwell coupling
     :exclude-members: options, allocate
     :show-inheritance:
 
+.. automodule:: struphy.propagators.vlasov_maxwell_coupling
+    :members:
+    :exclude-members: options, allocate
+    :show-inheritance:
+
 MHD-kinetic coupling
 ^^^^^^^^^^^^^^^^^^^^
 

examples/VlasovMaxwellQuasiNeutral/weibel_quasi_neutral/params_weibel_quasi_neutral.py

@@ -0,0 +1,195 @@
+"""Weibel instability with the VlasovMaxwellQuasiNeutral model, with and without displacement current
+(`--lam 1` or `--lam 0`), at a time step 20 times larger than in the VlasovMaxwellOneSpecies example."""
+
+import argparse
+from pathlib import Path
+
+# Command line options
+# no -h here when imported by the pproc script, which has its own help
+parser = argparse.ArgumentParser(add_help=__name__ == "__main__")
+parser.add_argument(
+    "--out-folder",
+    default=str(Path(__file__).resolve().parent),
+    help="Folder in which the output folder is created (default: beside this file).",
+)
+parser.add_argument("--sim-folder", default=None, help="Name of the output folder (default: sim_lam<lam>).")
+parser.add_argument(
+    "--lam",
+    type=float,
+    default=0.0,
+    help="lam = Omega_c x / c, lam^2 multiplies the displacement current (0: quasi-neutral).",
+)
+parser.add_argument("--dt", type=float, default=1.0, help="Time step.")
+parser.add_argument("--tend", type=float, default=400.0, help="End time.")
+parser.add_argument(
+    "--test-run",
+    action="store_true",
+    help="Few markers and three saved time steps, for the tests (test_examples.py).",
+)
+args, _ = parser.parse_known_args()
+if args.sim_folder is None:
+    args.sim_folder = f"sim_lam{args.lam:g}"
+if args.test_run:
+    args.tend = 2 * args.dt
+
+name = "Weibel instability in the quasi-neutral limit"
+
+description = f"""
+Weibel instability with the VlasovMaxwellQuasiNeutral model: 1D-3V kinetic electrons (PIC) with Maxwell's
+equations in a periodic box with a neutralizing ion background, solved with the implicit, energy conserving
+VlasovMaxwellCoupling step. The electrons are hotter in v2 than in v1 (along the wave vector); this temperature
+anisotropy makes a small transverse magnetic field B3 grow exponentially until magnetic trapping saturates it.
+Here lam = {args.lam:g}.
+
+Physics. The model uses the cyclotron normalization (time in units of 1/Omega_c, velocities in units of Omega_c x),
+with the equation parameters lam = Omega_c x / c and kappa = Omega_p x / c. With kappa = 1 lengths are in units of
+the skin depth c/omega_p; with lam = 1 also Omega_c = omega_p, so time is in units of 1/omega_p and velocities in
+units of c, as in the VlasovMaxwellOneSpecies Weibel example (alpha = 1). With epsilon = -1 (electrons) the model
+solves
+    df/dt + v . grad_x f - (E + v x B) . grad_v f = 0,
+    lam^2 dE/dt = curl B + kappa^2 int v f d^3v,
+    dB/dt = -curl E.
+lam = 1 is the Vlasov-Maxwell system; lam = 0 drops the displacement current (quasi-neutral, Darwin-type limit),
+where E is the Lagrange multiplier of the constraint curl B = J. For a bi-Maxwellian with thermal speeds vth1 (along
+k) and vth2, the transverse mode (E2, B3) obeys (kappa = 1)
+    lam^2 omega^2 - k^2 - 1 + (vth2^2 / vth1^2) [1 + xi Z(xi)] = 0,   xi = omega / (sqrt(2) k vth1)
+(the 1 is the plasma term, which stays for lam = 0). The mode is purely growing, omega = i gamma, and
+lam^2 gamma^2 ~ 8e-4 is negligible against k^2 + 1 = 2.56: the Weibel instability is a low-frequency, transverse
+(div J = 0) mode that the quasi-neutral limit keeps.
+
+Setup. As in the VlasovMaxwellOneSpecies Weibel example: periodic box of length 2 pi / k with k = 1.25, 32 cells
+and cubic splines in x, Maxwellian electrons with n = 1 and vth1 = 0.02 / sqrt(2), vth2 = sqrt(12) vth1 (v3 = 0),
+the seed field B3 = -1e-4 cos(k x), 100000 markers with full-f weights (the quasi-neutral limit needs full-f
+weights), loaded by an antithetic Sobol sequence (each Sobol point reflected in all phase-space directions, so that
+the mean velocities of each group of markers vanish), which removes most of the noise of the current. Time step
+dt = {args.dt:g} (20 times the dt = 0.05 of the VlasovMaxwellOneSpecies example; the implicit coupling step is stable
+for dt omega_p > 1) up to t = {args.tend:g}, Lie-Trotter splitting.
+
+Expected result. For k = 1.25 the dispersion relation gives gamma = 0.02784 for both lam = 1 and lam = 0 (the
+growth rate of the amplitude of B3; the magnetic energy grows as exp(2 gamma t)). The total energy en_tot conserved by the
+coupling step changes only by the splitting with the particle pushers. The pproc fits gamma to the magnetic energy
+(between t = Tend/5 and 2 Tend/5), compares with the analytical value and plots the energies. Measured with this
+setup (serial, 2026-10-23): gamma = 0.02805 for lam = 1 and gamma = 0.02758 for lam = 0, within 1 % of the
+analytical value; the relative change of en_tot stays below 6e-9. With pseudo-random markers instead, the fitted rate
+is about 10 % low and scatters by about 7 % between seeds (marker noise); the large time step changes it only by
+about 1 %. The verification test test_verif_VlasovMaxwellQuasiNeutral.py runs this setup for lam = 0 up to
+t = 160.
+
+For lam = 0 the coupling step projects the initial state onto curl B = J (option initial_projection, default):
+the constraint holds only at the mid-point of each step, so initial data that violate it (here by the marker noise
+of J) would oscillate with period 2.
+
+References:
+- E. S. Weibel, "Spontaneously growing transverse waves in a plasma due to an anisotropic velocity distribution",
+  Phys. Rev. Lett. 2, 83 (1959).
+- M. Kraus, K. Kormann, P. J. Morrison and E. Sonnendruecker, "GEMPIC: geometric electromagnetic particle-in-cell
+  methods", J. Plasma Phys. 83, 905830401 (2017).
+"""
+
+# ------------------
+# Import Struphy API
+# ------------------
+
+import logging
+import math
+
+from struphy import (
+    BaseUnits,
+    BoundaryParameters,
+    DerhamOptions,
+    EnvironmentOptions,
+    LoadingParameters,
+    SavingParameters,
+    Simulation,
+    SortingParameters,
+    Time,
+    WeightsParameters,
+    domains,
+    grids,
+    maxwellians,
+    perturbations,
+    set_logging_level,
+)
+
+set_logging_level(logging.WARNING)
+
+# ---------------------
+# Instance of the model
+# ---------------------
+
+from struphy.models import VlasovMaxwellQuasiNeutral
+
+model = VlasovMaxwellQuasiNeutral(base_units=BaseUnits(), epsilon=-1.0, lam=args.lam, kappa=1.0)
+
+k = 1.25
+B_pert_amp = -1e-4
+vth1_background_val = 0.02 / math.sqrt(2)
+vth2_background_val = vth1_background_val * math.sqrt(12)
+
+model.em_fields.e_field.save_data = True
+model.em_fields.b_field.save_data = True
+model.kinetic_ions.var.save_data = False
+
+# --------------------------
+# Instance of the simulation
+# --------------------------
+
+env = EnvironmentOptions(out_folders=args.out_folder, sim_folder=args.sim_folder)
+time_opts = Time(dt=args.dt, Tend=args.tend, split_algo="LieTrotter")
+Lx = 2 * math.pi / k
+domain = domains.Cuboid(r1=Lx)
+grid = grids.TensorProductGrid(num_elements=(32, 1, 1))
+derham_opts = DerhamOptions(degree=(3, 1, 1))
+
+sim = Simulation(
+    model=model,
+    name=name,
+    description=description,
+    params_path=__file__,
+    env=env,
+    time_opts=time_opts,
+    domain=domain,
+    grid=grid,
+    derham_opts=derham_opts,
+)
+
+# -------------------
+# Particle parameters
+# -------------------
+
+model.kinetic_ions.set_markers(
+    loading_params=LoadingParameters(
+        Np=2000 if args.test_run else 100000,
+        # quasi-random markers reflected in all phase-space directions (no seed): much less noise in the current
+        loading="sobol_antithetic",
+        set_zero_velocity=(False, False, True),
+        moments=(0.0, 0.0, 0.0, vth1_background_val, vth2_background_val, 1.0),
+    ),
+    weights_params=WeightsParameters(control_variate=False),
+    boundary_params=BoundaryParameters(),
+    sorting_params=SortingParameters(boxes_per_dim=(16, 1, 1), do_sort=True),
+    saving_params=SavingParameters(n_markers=0),
+    bufsize=2.0,
+)
+
+# ------------------
+# Propagator options
+# ------------------
+
+model.propagators.coupling.options = model.propagators.coupling.Options()
+model.propagators.push_eta.options = model.propagators.push_eta.Options()
+model.propagators.push_vxb.options = model.propagators.push_vxb.Options()
+model.initial_poisson.options = model.initial_poisson.Options(stab_mat="M0")
+
+# ------------------
+# Initial conditions
+# ------------------
+
+model.kinetic_ions.var.add_background(
+    maxwellians.Maxwellian3D(vth1=(vth1_background_val, None), vth2=(vth2_background_val, None))
+)
+model.em_fields.b_field.add_perturbation(perturbation=perturbations.ModesCos(amps=(B_pert_amp,), ls=(1,), comp=2))
+
+if __name__ == "__main__":
+    # post-processing and plots: pproc_weibel_quasi_neutral.py
+    sim.run()

examples/VlasovMaxwellQuasiNeutral/weibel_quasi_neutral/pproc_weibel_quasi_neutral.py

@@ -0,0 +1,71 @@
+"""Post-process the quasi-neutral Weibel example: growth rate of the magnetic energy and energy conservation.
+
+Run as ``python pproc_weibel_quasi_neutral.py [--lam L]`` (same flags as the params file) or with the path of the
+output folder. Writes ``energies.png`` into ``<sim_folder>/results`` and prints the fitted growth rate.
+"""
+
+import argparse
+import sys
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from struphy import Output
+
+SCRIPT_DIR = Path(__file__).resolve().parent
+sys.path.insert(0, str(SCRIPT_DIR))
+from params_weibel_quasi_neutral import args  # noqa: E402
+
+DEFAULT_OUTPUT = Path(args.out_folder) / args.sim_folder
+GAMMA_EXACT = 0.02784
+
+
+def main(path_out=DEFAULT_OUTPUT):
+    out = Output(path_out)
+    scalars = out.scalars
+    t = scalars["en_B"].t.values
+    en_B = scalars["en_B"].values
+    en_tot = scalars["en_tot"].values
+
+    results_dir = Path(path_out) / "results"
+    results_dir.mkdir(exist_ok=True)
+
+    # growth rate from the linear phase, Tend/5 < t < 2 Tend/5: en_B ~ exp(2 gamma t)
+    mask = (t > t[-1] / 5) & (t < 2 * t[-1] / 5)
+    if np.count_nonzero(mask) > 2:
+        slope, offset = np.polyfit(t[mask], np.log(en_B[mask]), 1)
+        gamma = slope / 2
+        print(
+            f"fitted gamma = {gamma:.5f}, analytical gamma = {GAMMA_EXACT} (rel. error {gamma / GAMMA_EXACT - 1:.2%})"
+        )
+    else:
+        gamma = None
+        print("Too few saved time steps for a growth-rate fit.")
+    print(f"max relative change of en_tot: {np.max(np.abs(en_tot - en_tot[0])) / abs(en_tot[0]):.2e}")
+
+    fig, axs = plt.subplots(1, 2, figsize=(11, 4))
+    axs[0].semilogy(t, en_B, label="en_B")
+    axs[0].semilogy(t, scalars["en_E"].values, label="en_E")
+    if gamma is not None:
+        axs[0].semilogy(t[mask], np.exp(offset + slope * t[mask]), "k--", label=f"fit, gamma = {gamma:.4f}")
+        axs[0].semilogy(
+            t, en_B[mask][0] * np.exp(2 * GAMMA_EXACT * (t - t[mask][0])), "r:", label=f"gamma = {GAMMA_EXACT}"
+        )
+    axs[0].set_ylim(en_B.min() / 10, en_B.max() * 10)
+    axs[0].set_xlabel("t")
+    axs[0].legend()
+    axs[1].plot(t, (en_tot - en_tot[0]) / abs(en_tot[0]))
+    axs[1].set_xlabel("t")
+    axs[1].set_ylabel("(en_tot - en_tot(0)) / en_tot(0)")
+    fig.suptitle(f"Weibel instability, {out.path_out.name}")
+    fig.tight_layout()
+    fig.savefig(results_dir / "energies.png")
+    plt.close(fig)
+
+
+if __name__ == "__main__":
+    pproc_parser = argparse.ArgumentParser(description="Post-process a saved simulation run.")
+    pproc_parser.add_argument("path_out", nargs="?", default=DEFAULT_OUTPUT, help="Output folder of the run.")
+    pproc_args, _ = pproc_parser.parse_known_args()
+    main(pproc_args.path_out)

src/struphy/api/models/__init__.py

@@ -46,6 +46,7 @@
 from struphy.models.vlasov import Vlasov
 from struphy.models.vlasov_ampere_one_species import VlasovAmpereOneSpecies
 from struphy.models.vlasov_maxwell_one_species import VlasovMaxwellOneSpecies
+from struphy.models.vlasov_maxwell_quasi_neutral import VlasovMaxwellQuasiNeutral
 from struphy.models.vlasov_poisson_adiabatic import VlasovPoissonAdiabatic
 from struphy.models.vlasov_poisson_adiabatic_pauli import VlasovPoissonAdiabaticPauli
 
@@ -90,6 +91,7 @@
     "Vlasov",
     "VlasovAmpereOneSpecies",
     "VlasovMaxwellOneSpecies",
+    "VlasovMaxwellQuasiNeutral",
     "VlasovPoissonAdiabatic",
     "VlasovPoissonAdiabaticPauli",
 ]

src/struphy/models/AGENTS.md

@@ -23,7 +23,7 @@ Read this before creating or renaming a model in `src/struphy/models/`. The gene
   0. `self.params = copy.deepcopy(locals())`. All constructor arguments need defaults, and only
      bool/int/float/str/tuple/list/None/`BaseUnits` values survive serialization.
   1. instantiate the species as attributes;
-  2. `self.setup_equation_params(base_units=base_units)` (derives units and `alpha`, `epsilon`, `kappa`);
+  2. `self.setup_equation_params(base_units=base_units)` (derives units and `alpha`, `epsilon`, `lam`, `kappa`);
   3. `self.propagators = self.Propagators()`;
   4. link variables: `self.propagators.<prop>.variables.<var> = self.<species>.<var>`;
   5. `self.scalars = Scalars(...)` with the tracked quantities (`scalars.py`, e.g. `BilinearEnergyFEEC`,
@@ -44,6 +44,8 @@ Read this before creating or renaming a model in `src/struphy/models/`. The gene
 - All 8 doc classmethods, written as reST docstrings with math in `.. math::` / `:math:`:
   `doc_pde`, `doc_normalization`, `doc_scalar_quantities`, `doc_long_description`, `doc_examples`,
   `doc_use_cases`, `doc_cannot_be_used_for`, and `doc_discretization`.
+  Put `doc_normalization` directly below `doc_pde` (also in subclasses that change the normalization), since the
+  normalization defines the parameters appearing in the PDEs.
   `doc_discretization` must **return** a string: list the propagators in splitting order with `:class:` links and
   embed each propagator's `__doc__` (see `maxwell.py`). `test_all_models_expose_new_doc_api` fails if any is missing or empty.
 - Add an `.. autoclass:: struphy.models.<Name>` block (same options as its neighbours) under the matching category in

src/struphy/models/__init__.py

@@ -49,6 +49,7 @@
     "Vlasov": "struphy.models.vlasov",
     "VlasovAmpereOneSpecies": "struphy.models.vlasov_ampere_one_species",
     "VlasovMaxwellOneSpecies": "struphy.models.vlasov_maxwell_one_species",
+    "VlasovMaxwellQuasiNeutral": "struphy.models.vlasov_maxwell_quasi_neutral",
     "VlasovPoissonAdiabatic": "struphy.models.vlasov_poisson_adiabatic",
     "VlasovPoissonAdiabaticPauli": "struphy.models.vlasov_poisson_adiabatic_pauli",
 }
@@ -94,6 +95,7 @@
     from struphy.models.vlasov import Vlasov
     from struphy.models.vlasov_ampere_one_species import VlasovAmpereOneSpecies
     from struphy.models.vlasov_maxwell_one_species import VlasovMaxwellOneSpecies
+    from struphy.models.vlasov_maxwell_quasi_neutral import VlasovMaxwellQuasiNeutral
     from struphy.models.vlasov_poisson_adiabatic import VlasovPoissonAdiabatic
     from struphy.models.vlasov_poisson_adiabatic_pauli import VlasovPoissonAdiabaticPauli
 
@@ -138,6 +140,7 @@
     "Vlasov",
     "VlasovAmpereOneSpecies",
     "VlasovMaxwellOneSpecies",
+    "VlasovMaxwellQuasiNeutral",
     "VlasovPoissonAdiabatic",
     "VlasovPoissonAdiabaticPauli",
 ]