diff --git a/doc/sections/subsections/models-all.rst b/doc/sections/subsections/models-all.rst index 9b5ad9a5f..b896f7082 100644 --- a/doc/sections/subsections/models-all.rst +++ b/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: diff --git a/doc/sections/subsections/models-normalization.rst b/doc/sections/subsections/models-normalization.rst index 058d5d801..a204334bf 100644 --- a/doc/sections/subsections/models-normalization.rst +++ b/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. diff --git a/doc/sections/subsections/propagators-avail.rst b/doc/sections/subsections/propagators-avail.rst index 7fd276e1c..ac971567a 100644 --- a/doc/sections/subsections/propagators-avail.rst +++ b/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 ^^^^^^^^^^^^^^^^^^^^ diff --git a/examples/VlasovMaxwellQuasiNeutral/weibel_quasi_neutral/params_weibel_quasi_neutral.py b/examples/VlasovMaxwellQuasiNeutral/weibel_quasi_neutral/params_weibel_quasi_neutral.py new file mode 100644 index 000000000..5b339fb28 --- /dev/null +++ b/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).") +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() diff --git a/examples/VlasovMaxwellQuasiNeutral/weibel_quasi_neutral/pproc_weibel_quasi_neutral.py b/examples/VlasovMaxwellQuasiNeutral/weibel_quasi_neutral/pproc_weibel_quasi_neutral.py new file mode 100644 index 000000000..aff17997c --- /dev/null +++ b/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 ``/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) diff --git a/src/struphy/api/models/__init__.py b/src/struphy/api/models/__init__.py index d77c89bf6..3a390184c 100644 --- a/src/struphy/api/models/__init__.py +++ b/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", ] diff --git a/src/struphy/models/AGENTS.md b/src/struphy/models/AGENTS.md index cb7a15b22..3ca957989 100644 --- a/src/struphy/models/AGENTS.md +++ b/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..variables. = self..`; 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.` block (same options as its neighbours) under the matching category in diff --git a/src/struphy/models/__init__.py b/src/struphy/models/__init__.py index e3fd8abbc..ba12868c5 100644 --- a/src/struphy/models/__init__.py +++ b/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", ] diff --git a/src/struphy/models/species.py b/src/struphy/models/species.py index 3c8b47657..b43464713 100644 --- a/src/struphy/models/species.py +++ b/src/struphy/models/species.py @@ -56,11 +56,14 @@ class Species(metaclass=ABCMeta): epsilon : float, optional Normalized cyclotron period: 1 / (cyclotron frequency × time unit). If None, computed from units and charge/mass numbers. + lam : float, optional + Cyclotron frequency × length unit / speed of light. + If None, computed from units and charge/mass numbers. kappa : float, optional - Normalized plasma frequency: plasma frequency × time unit. + Plasma frequency × length unit / speed of light. If None, computed from units and charge/mass numbers. equation_params : EquationParameters - Object containing normalization parameters (alpha, epsilon, kappa) for scaled equations. + Object containing normalization parameters (alpha, epsilon, lam, kappa) for scaled equations. Methods ------- @@ -94,6 +97,7 @@ def to_dict(self) -> dict: "mass_number": self.mass_number, "alpha": self.alpha, "epsilon": self.epsilon, + "lam": self.lam, "kappa": self.kappa, "variables": {name: variable.to_dict() for name, variable in self.variables.items()}, } @@ -105,11 +109,12 @@ def init_variables( alpha: float = None, epsilon: float = None, kappa: float = None, + lam: float = None, ): """Create variables dict and set physical and equation parameters for a plasma species. Sets the charge and mass numbers, and optionally overrides normalized equation parameters - (alpha, epsilon, kappa) that would otherwise be computed from physical units. + (alpha, epsilon, lam, kappa) that would otherwise be computed from physical units. Parameters ---------- @@ -124,7 +129,10 @@ def init_variables( Normalized cyclotron period: 1 / (cyclotron frequency × time unit). If None, computed from units and charge/mass numbers (default = None). kappa : float, optional - Normalized plasma frequency: plasma frequency × time unit. + Plasma frequency × length unit / speed of light. + If None, computed from units and charge/mass numbers (default = None). + lam : float, optional + Cyclotron frequency × length unit / speed of light. If None, computed from units and charge/mass numbers (default = None).""" # create variables dict @@ -141,6 +149,7 @@ def init_variables( self._alpha = alpha self._epsilon = epsilon self._kappa = kappa + self._lam = lam @property def variables(self) -> dict: @@ -174,9 +183,16 @@ def epsilon(self) -> float: self._epsilon = None return self._epsilon + @property + def lam(self) -> float: + """Cyclotron frequency times length unit over speed of light, Omega_c * x_unit / c (default = None).""" + if not hasattr(self, "_lam"): + self._lam = None + return self._lam + @property def kappa(self) -> float: - """The normalized plasma frequency, Omega_p * time_unit (default = None).""" + """Plasma frequency times length unit over speed of light, Omega_p * x_unit / c (default = None).""" if not hasattr(self, "_kappa"): self._kappa = None return self._kappa @@ -191,6 +207,7 @@ def __init__( alpha: float = None, epsilon: float = None, kappa: float = None, + lam: float = None, ): self.species = species @@ -221,8 +238,17 @@ def __init__( if MPI.COMM_WORLD.Get_rank() == 0: warnings.warn(f"Override equation parameter {self.epsilon =}") + # lam and kappa: the coefficients lam^2 of the displacement current and kappa^2 of the plasma current in + # Ampere's law for the cyclotron velocity scale (kappa = alpha * lam) + if lam is None: + self.lam = om_c * units.x / con.c + else: + self.lam = lam + if MPI.COMM_WORLD.Get_rank() == 0: + warnings.warn(f"Override equation parameter {self.lam =}") + if kappa is None: - self.kappa = om_p * units.t + self.kappa = om_p * units.x / con.c else: self.kappa = kappa if MPI.COMM_WORLD.Get_rank() == 0: @@ -243,7 +269,8 @@ def setup_equation_params(self, units: Units): * alpha = plasma-frequenca / cyclotron frequency * epsilon = 1 / (cyclotron frequency * time unit) - * kappa = plasma frequency * time unit + * lam = cyclotron frequency * length unit / speed of light + * kappa = plasma frequency * length unit / speed of light """ self._equation_params = self.EquationParameters( species=self, @@ -251,6 +278,7 @@ def setup_equation_params(self, units: Units): alpha=self.alpha, epsilon=self.epsilon, kappa=self.kappa, + lam=self.lam, ) diff --git a/src/struphy/models/tests/verification/test_verif_VlasovMaxwellQuasiNeutral.py b/src/struphy/models/tests/verification/test_verif_VlasovMaxwellQuasiNeutral.py new file mode 100644 index 000000000..c41e86b9c --- /dev/null +++ b/src/struphy/models/tests/verification/test_verif_VlasovMaxwellQuasiNeutral.py @@ -0,0 +1,116 @@ +import logging +import math +import os +import shutil + +import h5py +import numpy as np +from maybempi import MPI + +from struphy import ( + BoundaryParameters, + DerhamOptions, + EnvironmentOptions, + LoadingParameters, + SavingParameters, + Simulation, + SortingParameters, + Time, + WeightsParameters, + domains, + grids, + maxwellians, + perturbations, +) +from struphy.models import VlasovMaxwellQuasiNeutral + +logger = logging.getLogger("struphy") + +# growth rate of the amplitude of B3 for k = 1.25, vth2^2 / vth1^2 = 12 (Kraus et al. 2017); the same for lam = 0, +# since the displacement-current term lam^2 gamma^2 ~ 8e-4 of the lam = 1 dispersion relation is << k^2 + 1 +GAMMA_EXACT = 0.02784 + + +def _weibel_growth_rate( + lam: float, Np: int, Tend: float, t_fit: tuple, dt: float = 1.0, loading: str = "sobol_antithetic" +) -> float: + """Runs the Weibel instability of the example weibel_quasi_neutral (1D-3V, k = 1.25, full-f markers) and returns + the growth rate fitted to the magnetic energy, en_B ~ exp(2 gamma t), for t_fit[0] < t < t_fit[1]. + + The antithetic Sobol loading (each Sobol point reflected u -> 1 - u in all phase-space directions, 64 markers + per point, so that the mean velocities of each group vanish) removes most of the noise of the current, which + limits the fit: with pseudo-random markers the fitted rate scatters by about 7 % between seeds + (1e5 markers).""" + # electrons, cyclotron normalization with Omega_c = Omega_p (kappa = 1): lengths in units of the skin depth + model = VlasovMaxwellQuasiNeutral(epsilon=-1.0, lam=lam, kappa=1.0) + + test_folder = os.path.join(os.getcwd(), "struphy_verification_tests") + out_folders = os.path.join(test_folder, "VlasovMaxwellQuasiNeutral") + env = EnvironmentOptions(out_folders=out_folders, sim_folder=f"weibel_lam{lam:g}_Np{Np}_{loading}") + + k = 1.25 + vth1 = 0.02 / math.sqrt(2) + vth2 = vth1 * math.sqrt(12) + + model.kinetic_ions.set_markers( + loading_params=LoadingParameters( + Np=Np, + loading=loading, + set_zero_velocity=(False, False, True), + moments=(0.0, 0.0, 0.0, vth1, vth2, 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, + ) + 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") + + model.kinetic_ions.var.add_background(maxwellians.Maxwellian3D(vth1=(vth1, None), vth2=(vth2, None))) + model.em_fields.b_field.add_perturbation(perturbation=perturbations.ModesCos(amps=(-1e-4,), ls=(1,), comp=2)) + + sim = Simulation( + model=model, + env=env, + time_opts=Time(dt=dt, Tend=Tend, split_algo="LieTrotter"), + domain=domains.Cuboid(r1=2 * math.pi / k), + grid=grids.TensorProductGrid(num_elements=(32, 1, 1)), + derham_opts=DerhamOptions(degree=(3, 1, 1)), + ) + sim.run() + + # the scalars are global (reduced over the ranks), so the file of rank 0 holds them for all ranks + comm = MPI.COMM_WORLD + comm.Barrier() + with h5py.File(os.path.join(env.path_out, "data", "data_proc0.hdf5"), "r") as file: + t = np.asarray(file["time/value"]) + en_B = np.asarray(file["scalar/en_B"]) + comm.Barrier() + + mask = (t > t_fit[0]) & (t < t_fit[1]) + gamma = np.polyfit(t[mask], np.log(en_B[mask]), 1)[0] / 2 + return float(gamma) + + +def test_weibel_growth_rate_quasi_neutral(): + """Weibel instability in the quasi-neutral limit (lam = 0, no displacement current): the magnetic energy grows + with the analytical growth rate of the Vlasov-Maxwell system, gamma = 0.02784, which the displacement current does + not change. With 1e5 full-f markers (antithetic Sobol loading) and dt = 1 the fitted rate is close to the + analytical one (about 0.5 %), on 1, 2 and 4 MPI ranks.""" + gamma = _weibel_growth_rate(lam=0.0, Np=100_000, Tend=160.0, t_fit=(80.0, 150.0)) + logger.info(f"Weibel growth rate (lam = 0): {gamma = }, analytical {GAMMA_EXACT}") + assert abs(gamma / GAMMA_EXACT - 1) < 0.07, f"{gamma = } vs. {GAMMA_EXACT = }" + + comm = MPI.COMM_WORLD + comm.Barrier() + if comm.Get_rank() == 0: + shutil.rmtree(os.path.join(os.getcwd(), "struphy_verification_tests"), ignore_errors=True) + comm.Barrier() + + +if __name__ == "__main__": + test_weibel_growth_rate_quasi_neutral() diff --git a/src/struphy/models/vlasov_maxwell_one_species.py b/src/struphy/models/vlasov_maxwell_one_species.py index 66ae1185d..46573cc61 100644 --- a/src/struphy/models/vlasov_maxwell_one_species.py +++ b/src/struphy/models/vlasov_maxwell_one_species.py @@ -116,12 +116,7 @@ def __init__( self.propagators = self.Propagators(b2_var=self.em_fields.b_field) # 4. assign variables to propagators - self.propagators.maxwell.variables.e = self.em_fields.e_field - self.propagators.maxwell.variables.b = self.em_fields.b_field - self.propagators.push_eta.variables.var = self.kinetic_ions.var - self.propagators.push_vxb.variables.ions = self.kinetic_ions.var - self.propagators.coupling_va.variables.e = self.em_fields.e_field - self.propagators.coupling_va.variables.ions = self.kinetic_ions.var + self._link_propagator_variables() # 5. define scalars to be tracked during simulation electric_energy = BilinearEnergyFEEC(self.em_fields.e_field) @@ -160,6 +155,18 @@ def __init__( # property to measure violation of gauss law from control variate self.measure_gauss_law = measure_gauss_law + def _link_propagator_variables(self): + """Assign the species variables to the propagators (overridden by models with other propagators). + + :meta private: + """ + self.propagators.maxwell.variables.e = self.em_fields.e_field + self.propagators.maxwell.variables.b = self.em_fields.b_field + self.propagators.push_eta.variables.var = self.kinetic_ions.var + self.propagators.push_vxb.variables.ions = self.kinetic_ions.var + self.propagators.coupling_va.variables.e = self.em_fields.e_field + self.propagators.coupling_va.variables.ions = self.kinetic_ions.var + @property def bulk_species(self): return self.kinetic_ions @@ -251,17 +258,20 @@ def calculate_gauss_error(self): ## default parameters def generate_default_parameter_file(self, path=None, prompt=True): params_path = super().generate_default_parameter_file(path=path, prompt=prompt) + # the options of the initial Poisson solve go after those of the last propagator + last_prop = list(self.propagators.__dict__)[-1] + model_name = type(self).__name__ new_file = [] with open(params_path, "r") as f: for line in f: - if "coupling_va.Options" in line: + if f"{last_prop}.Options" in line: new_file += [line] new_file += ["model.initial_poisson.options = model.initial_poisson.Options()\n"] elif "saving_params = " in line: new_file += ["\nbinplot = BinningPlot(slice='e1', n_bins=128, ranges=(0.0, 1.0))\n"] new_file += ["saving_params = SavingParameters(binning_plots=(binplot,))\n\n"] - elif "VlasovMaxwellOneSpecies()" in line: - new_file += ["\nmodel = VlasovMaxwellOneSpecies(measure_gauss_law=True)\n"] + elif f"{model_name}()" in line: + new_file += [f"\nmodel = {model_name}(measure_gauss_law=True)\n"] else: new_file += [line] diff --git a/src/struphy/models/vlasov_maxwell_quasi_neutral.py b/src/struphy/models/vlasov_maxwell_quasi_neutral.py new file mode 100644 index 000000000..95b0a0341 --- /dev/null +++ b/src/struphy/models/vlasov_maxwell_quasi_neutral.py @@ -0,0 +1,328 @@ +import copy +import logging + +from struphy import BaseUnits +from struphy.models.scalars import FunctionScalarFEEC, Scalars +from struphy.models.species import ParticleSpecies +from struphy.models.variables import PICVariable +from struphy.models.vlasov_maxwell_one_species import VlasovMaxwellOneSpecies +from struphy.pic.accumulation.kernels.charge_density_0form import charge_density_0form +from struphy.pic.accumulation.particles_to_grid import ParticlesToGrid +from struphy.propagators.poisson_solve import PoissonSolve +from struphy.propagators.push_eta import PushEta +from struphy.propagators.push_vxb import PushVxB +from struphy.propagators.vlasov_maxwell_coupling import VlasovMaxwellCoupling + +logger = logging.getLogger("struphy") + + +class VlasovMaxwellQuasiNeutral(VlasovMaxwellOneSpecies): + """Vlasov-Maxwell equations for one kinetic species in the cyclotron normalization, with an implicit, energy + conserving field-particle coupling that also holds in the quasi-neutral limit (no displacement current). + + Same equations as :class:`~struphy.models.vlasov_maxwell_one_species.VlasovMaxwellOneSpecies`, normalized with the + cyclotron frequency (velocity scale ``"cyclotron"``), so that the coefficients :math:`\\lambda^2` of the + displacement current and :math:`\\kappa^2` of the plasma current appear in Ampère's law. Ampère's law, Faraday's + law and the electric force on the particles are solved together in one Crank-Nicolson step + (:class:`~struphy.propagators.vlasov_maxwell_coupling.VlasovMaxwellCoupling`); ``lam=0.0`` is the quasi-neutral + limit. + + Parameters + ---------- + base_units: BaseUnits + Base units for normalization (default: BaseUnits()) + charge_number: int + Charge number (in units of the positive elementary charge) of the species (default: 1) + mass_number: float + Mass number (in units of Proton mass) of the species (default: 1.0) + alpha: float, optional + Dimensionless parameter: plasma frequency / cyclotron frequency. If None, computed from units and charge/mass numbers. + It does not enter the equations of this model (which use lam and kappa = alpha lam); when overriding lam and + kappa, alpha is not adapted. + epsilon: float, optional + Normalized cyclotron period: 1 / (cyclotron frequency × time unit). If None, computed from units and charge/mass + numbers (equal to 1 in this model's normalization). + lam: float, optional + Cyclotron frequency × length unit / speed of light (coefficient lam^2 of the displacement current; 0 is the + quasi-neutral limit). If None, computed from units and charge/mass numbers. + kappa: float, optional + Plasma frequency × length unit / speed of light (coefficient kappa^2 of the plasma current). If None, computed + from units and charge/mass numbers. + measure_gauss_law: bool + Whether to track the Gauss-law error as a scalar quantity (default: False; not available for lam = 0) + """ + + ## species + + class KineticIons(ParticleSpecies): + def __init__( + self, + charge_number: int = 1, + mass_number: float = 1.0, + alpha: float = None, + epsilon: float = None, + lam: float = None, + kappa: float = None, + ): + self.var = PICVariable(space="Particles6D") + self.init_variables( + charge_number=charge_number, + mass_number=mass_number, + alpha=alpha, + epsilon=epsilon, + lam=lam, + kappa=kappa, + ) + + ## propagators + + class Propagators: + def __init__( + self, + b2_var=None, + ): + self.coupling = VlasovMaxwellCoupling() + self.push_eta = PushEta() + self.push_vxb = PushVxB(b2_var=b2_var) + + ## abstract methods + + def __init__( + self, + base_units: BaseUnits = BaseUnits(), + charge_number: int = 1, + mass_number: float = 1.0, + alpha: float = None, + epsilon: float = None, + lam: float = None, + kappa: float = None, + measure_gauss_law: bool = False, + ): + # 0. store input parameters + self.params = copy.deepcopy(locals()) + + # 1. instantiate all species + self.em_fields = self.EMFields() + self.kinetic_ions = self.KineticIons(charge_number, mass_number, alpha, epsilon, lam, kappa) + + # 2. derive units (must be done after instantiating species to access charge and mass numbers) + self.setup_equation_params(base_units=base_units) + + # 3. instantiate all propagators + self.propagators = self.Propagators(b2_var=self.em_fields.b_field) + + # 4. assign variables to propagators + self._link_propagator_variables() + + # 5. define scalars to be tracked during simulation: the parts of the energy conserved by the coupling step + # (with the coefficients lam^2, c1, c2, c3 of the propagator, known after its allocation) + coupling = self.propagators.coupling + electric_energy = FunctionScalarFEEC(lambda: coupling.energies()[0]) + magnetic_energy = FunctionScalarFEEC(lambda: coupling.energies()[1]) + particle_energy = FunctionScalarFEEC(lambda: coupling.energies()[2]) + scalars_dict = { + "en_E": electric_energy, + "en_B": magnetic_energy, + "en_f": particle_energy, + "en_tot": electric_energy + magnetic_energy + particle_energy, + } + if measure_gauss_law: + # the MPI reduction (max over ranks) is done in calculate_gauss_error, not by summing the scalar + scalars_dict["gauss_error"] = FunctionScalarFEEC(self.calculate_gauss_error) + self.scalars = Scalars(**scalars_dict) + + # initial Poisson (not a propagator used in time stepping); its charge coefficient c2 / lam^2 is set in + # post_allocate, from the coefficients of the coupling propagator + particles_to_grid = ParticlesToGrid( + self.kinetic_ions.var, + "H1", + charge_density_0form, + ) + self.initial_poisson = PoissonSolve( + rho=particles_to_grid, + rho_coeffs=None, + delta_f_source=True, + ) + self.initial_poisson.variables.phi = self.em_fields.phi + + # property to measure violation of gauss law from control variate + self.measure_gauss_law = measure_gauss_law + + @property + def velocity_scale(self): + return "cyclotron" + + def _link_propagator_variables(self): + """Assign the species variables to the propagators. + + :meta private: + """ + self.propagators.coupling.variables.e = self.em_fields.e_field + self.propagators.coupling.variables.b = self.em_fields.b_field + self.propagators.coupling.variables.ions = self.kinetic_ions.var + self.propagators.push_eta.variables.var = self.kinetic_ions.var + self.propagators.push_vxb.variables.ions = self.kinetic_ions.var + + def post_allocate(self): + """Solve the initial Poisson equation lam^2 div E = c2 rho (Gauss' law consistent with the coupling step), + except in the quasi-neutral limit lam = 0, where the initial electric field (zero) does not enter the steps. + + :meta private: + """ + coupling = self.propagators.coupling + self._quasi_neutral = coupling._lam2 == 0.0 + if self._quasi_neutral: + logger.info( + "\nlam = 0 (quasi-neutral limit): no initial Poisson solve, the Gauss-law error is not measured." + ) + if self.measure_gauss_law: + logger.warning("measure_gauss_law: Gauss' law is not defined for lam = 0; gauss_error is NaN.") + return + + self.initial_poisson.rho_coeffs = coupling._c2 / coupling._lam2 + super().post_allocate() + + def calculate_gauss_error(self): + r"""Maximum norm of the weak Gauss-law residual + :math:`\mathbb G^\top \mathbb M^1 \mathbf e + \frac{c_2}{\lambda^2} \boldsymbol \rho`, see + :meth:`~struphy.models.vlasov_maxwell_one_species.VlasovMaxwellOneSpecies.calculate_gauss_error` + (NaN for :math:`\lambda = 0`).""" + if self._quasi_neutral: + return float("nan") + return super().calculate_gauss_error() + + @classmethod + def doc_pde(cls): + r"""**PDEs solved by model:** + + Vlasov equation: + + .. math:: + + \frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla f + \left( \mathbf{E} + \mathbf{v} \times \left( \mathbf{B} + \mathbf{B}_0 \right) \right) \cdot \frac{\partial f}{\partial \mathbf{v}} = 0 + + Ampère's law: + + .. math:: + + -\lambda^2 \frac{\partial \mathbf{E}}{\partial t} + \nabla \times \mathbf{B} = \kappa^2 \int_{\mathbb{R}^3} \mathbf{v} f \, \text{d}^3 \mathbf{v} + + Faraday's law: + + .. math:: + + \frac{\partial \mathbf{B}}{\partial t} + \nabla \times \mathbf{E} = 0 + + For :math:`\lambda = 0` (quasi-neutral limit) Ampère's law becomes :math:`\nabla \times \mathbf B = + \kappa^2 \int \mathbf v f \, \textrm d^3 \mathbf v`, which implies :math:`\nabla \cdot \mathbf J = 0`; the electric field is then determined by this + constraint, not by its own evolution. The coefficients of the coupling terms can be changed by the options + ``lam``, ``c1``, ``c2``, ``c3`` of the propagator + :class:`~struphy.propagators.vlasov_maxwell_coupling.VlasovMaxwellCoupling`. + + At initial time the weak Poisson equation (Gauss' law :math:`\lambda^2 \nabla \cdot \mathbf E = \kappa^2 \rho`) + is solved once for :math:`\lambda > 0`, + + .. math:: + + \lambda^2 \int_{\Omega} \nabla \psi^{\top} \cdot \nabla \phi \, \textrm{d} \mathbf{x} = \kappa^2 \int_{\Omega} \int_{\mathbb{R}^3} \psi \, (f - f_0) \, \text{d}^3 \mathbf{v} \, \textrm{d} \mathbf{x} \qquad \forall \ \psi \in H^1 + \\[2mm] + \mathbf{E}(t=0) = -\nabla \phi(t=0) + + (for :math:`\lambda = 0` the initial electric field does not enter the time stepping). Moreover, it is assumed + that :math:`\nabla \times \mathbf{B}_0 = \kappa^2 \int_{\mathbb{R}^3} \mathbf{v} f_0 \, \text{d}^3 \mathbf{v}`, + where :math:`\mathbf{B}_0` is the static equilibrium magnetic field. + """ + + @classmethod + def doc_normalization(cls): + r"""The model uses the cyclotron velocity of the species as reference velocity: + + .. math:: + + \hat v = \hat\Omega_c \hat x,\qquad \hat t = \frac{1}{\hat\Omega_c},\qquad \hat E = \hat B \hat v,\qquad \hat\phi = \hat E \hat x. + + The species parameters (see :ref:`equation_params`) are + + .. math:: + + \lambda = \frac{\hat\Omega_c \hat x}{c}\,,\qquad \kappa = \frac{\hat\Omega_p \hat x}{c} = \alpha \lambda\,, + + where :math:`c/\hat\Omega_c` and :math:`c/\hat\Omega_p` are the inertial lengths (:math:`\lambda` carries the + sign of the charge, like :math:`\hat\Omega_c`; only :math:`\lambda^2` enters). The Vlasov equation and + Faraday's law then have no coefficients, and Ampère's law + :math:`\frac{1}{c^2} \partial_t \mathbf E = \nabla \times \mathbf B - \mu_0 \mathbf J` becomes + :math:`\lambda^2 \partial_t \mathbf E = \nabla \times \mathbf B - \kappa^2 \mathbf J`. The quasi-neutral limit is + :math:`\lambda \to 0` at fixed :math:`\kappa`. With :math:`\lambda = \kappa = 1` the equations are those of :class:`~struphy.models.vlasov_maxwell_one_species.VlasovMaxwellOneSpecies` + with :math:`\alpha = 1`, i.e. lengths in units of :math:`c/\hat\Omega_p` and :math:`\hat\Omega_c = \hat\Omega_p`.""" + + @classmethod + def doc_scalar_quantities(cls): + r"""**The following scalars are tracked during simulation** (the parts of the energy conserved by the coupling + step, with the coefficients :math:`\lambda^2, c_1, c_2, c_3` of + :class:`~struphy.propagators.vlasov_maxwell_coupling.VlasovMaxwellCoupling`): + + - Electric field energy :math:`\frac{\lambda^2}{2} \int |\mathbf E|^2 \textrm d \mathbf x`: ``en_E`` + - Magnetic field energy :math:`\frac{c_1}{2} \int |\mathbf B|^2 \textrm d \mathbf x`: ``en_B`` + - Particle kinetic energy :math:`\frac{c_2}{2 c_3} \sum_p w_p |\mathbf v_p|^2`: ``en_f`` + (:math:`\frac{\kappa^2}{2} \sum_p w_p |\mathbf v_p|^2` for the default coefficients) + - Total energy: ``en_tot`` (conserved by the coupling step up to the solver tolerance) + - Optional Gauss-law diagnostic: ``gauss_error`` (NaN for :math:`\lambda = 0`)""" + + @classmethod + def doc_discretization(cls): + """Time integration is performed by the following propagators (in sequence): + + 1. :class:`~struphy.propagators.vlasov_maxwell_coupling.VlasovMaxwellCoupling` + 2. :class:`~struphy.propagators.push_eta.PushEta` + 3. :class:`~struphy.propagators.push_vxb.PushVxB` + """ + doc = rf"""**1. VlasovMaxwellCoupling:** + +{VlasovMaxwellCoupling.__doc__} + +**2. PushEta:** + +{PushEta.__doc__} + +**3. PushVxB:** + +{PushVxB.__doc__} +""" + return doc + + @classmethod + def doc_long_description(cls): + r"""VlasovMaxwellQuasiNeutral solves the Vlasov-Maxwell system in the cyclotron normalization, with one + implicit step for Maxwell's equations and the electric force, which conserves the energy exactly and allows time + steps larger than the inverse plasma frequency and the light-wave CFL limit. The coefficient :math:`\lambda^2` + of the displacement current can be set to zero (``lam=0``), which gives the quasi-neutral (Darwin-type) limit, + in which the electric field is the Lagrange multiplier of a divergence-free current.""" + + @classmethod + def doc_examples(cls): + r"""Create a quasi-neutral Vlasov-Maxwell model (lengths in units of the skin depth): + + .. code-block:: python + + from struphy.models import VlasovMaxwellQuasiNeutral + + model = VlasovMaxwellQuasiNeutral(lam=0.0, kappa=1.0) + model.em_fields.b_field + model.kinetic_ions.var + """ + + @classmethod + def doc_use_cases(cls): + r"""This model is appropriate for: + + - electromagnetic kinetic simulations with time steps above the plasma period + - the quasi-neutral limit of the Vlasov-Maxwell system (``lam=0.0``) + - comparison with :class:`~struphy.models.vlasov_maxwell_one_species.VlasovMaxwellOneSpecies`""" + + @classmethod + def doc_cannot_be_used_for(cls): + r"""This model is not suitable for: + + - the quasi-neutral limit with control-variate (delta-f) weights or with vacuum regions + - multi-species plasma dynamics without extension + - collisional kinetic closures""" diff --git a/src/struphy/propagators/__init__.py b/src/struphy/propagators/__init__.py index a4b34c03e..770e51e8b 100644 --- a/src/struphy/propagators/__init__.py +++ b/src/struphy/propagators/__init__.py @@ -50,6 +50,7 @@ "VariationalResistivity": "struphy.propagators.variational_resistivity", "VariationalViscosity": "struphy.propagators.variational_viscosity", "VlasovAmpereCoupling": "struphy.propagators.vlasov_ampere_coupling", + "VlasovMaxwellCoupling": "struphy.propagators.vlasov_maxwell_coupling", } if TYPE_CHECKING: # static analysis and IDEs see the eager imports @@ -102,6 +103,7 @@ from struphy.propagators.variational_resistivity import VariationalResistivity from struphy.propagators.variational_viscosity import VariationalViscosity from struphy.propagators.vlasov_ampere_coupling import VlasovAmpereCoupling + from struphy.propagators.vlasov_maxwell_coupling import VlasovMaxwellCoupling __all__ = [ "AdiabaticPhi", @@ -151,6 +153,7 @@ "VariationalResistivity", "VariationalViscosity", "VlasovAmpereCoupling", + "VlasovMaxwellCoupling", ] diff --git a/src/struphy/propagators/tests/test_matrix_free_6d.py b/src/struphy/propagators/tests/test_matrix_free_6d.py index 267c37362..51669e726 100644 --- a/src/struphy/propagators/tests/test_matrix_free_6d.py +++ b/src/struphy/propagators/tests/test_matrix_free_6d.py @@ -64,6 +64,48 @@ def run(matrix_free: bool): _compare(run) +@pytest.mark.parametrize("lam", [1.0, 0.0]) +def test_vlasov_maxwell_coupling_matrix_free(lam: float): + """VlasovMaxwellCoupling: e-field, b-field and velocities after two steps (in a curved domain), with and without + displacement current.""" + import numpy as np + + from struphy import DerhamOptions, LoadingParameters, Simulation, domains, grids, maxwellians, perturbations + from struphy.models import VlasovMaxwellQuasiNeutral + + def run(matrix_free: bool): + model = VlasovMaxwellQuasiNeutral(epsilon=-1.0, lam=lam, kappa=1.0) + model.kinetic_ions.set_markers(loading_params=LoadingParameters(Np=1000, seed=1234)) + model.kinetic_ions.var.add_background(maxwellians.Maxwellian3D()) + model.em_fields.b_field.add_perturbation( + perturbations.ModesSin(ls=(1,), ms=(1,), amps=(1e-2,), given_in_basis="2", comp=2) + ) + prop = model.propagators.coupling + prop.options = prop.Options(matrix_free=matrix_free, solver_params=_solver_params()) + Simulation( + model=model, + domain=domains.Colella(alpha=0.1), + grid=grids.TensorProductGrid(num_elements=(6, 4, 2)), + derham_opts=DerhamOptions(degree=(2, 2, 1)), + ).allocate() + particles = model.kinetic_ions.var.particles + # the reference density of the preconditioner, from the assembled or the matrix-free particle matrix + n_ref = np.array([prop._n_ref]) + # for lam=0 the change of b and V is the initial projection onto curl B = J; the steps then leave them + # unchanged and give e = e_mid = 0 (the initial e here), so e is not compared + changes = _steps( + prop, + lambda: ( + model.em_fields.e_field.spline.vector.toarray(), + model.em_fields.b_field.spline.vector.toarray(), + particles.markers[3:6, particles.valid_mks], + ), + ) + return [n_ref] + (changes if lam > 0.0 else changes[1:]) + + _compare(run) + + def test_efield_weights_coupling_matrix_free(): """EfieldWeightsCoupling: e-field and weights after two steps (in a curved domain).""" from struphy import DerhamOptions, LoadingParameters, Simulation, domains, grids, maxwellians, perturbations diff --git a/src/struphy/propagators/tests/test_vlasov_maxwell_coupling.py b/src/struphy/propagators/tests/test_vlasov_maxwell_coupling.py new file mode 100644 index 000000000..989b962cf --- /dev/null +++ b/src/struphy/propagators/tests/test_vlasov_maxwell_coupling.py @@ -0,0 +1,326 @@ +import logging + +import cunumpy as xp +import pytest + +from struphy import set_logging_level + +set_logging_level(logging.WARNING) + + +def _model( + model_name="VlasovMaxwellQuasiNeutral", control_variate=False, model_lam=1.0, measure_gauss_law=False, **options +): + """Allocated model (1D-ish slab, full-f electrons with alpha=1, epsilon=-1, a density and a B perturbation) and + its coupling propagator(s) with the given options.""" + from struphy import ( + DerhamOptions, + LoadingParameters, + Simulation, + WeightsParameters, + domains, + grids, + maxwellians, + models, + perturbations, + ) + from struphy.linear_algebra.solver import SolverParameters + + # electrons; lam = kappa = 1 in the cyclotron normalization equals alpha = 1 in the light normalization + if model_name == "VlasovMaxwellQuasiNeutral": + model = getattr(models, model_name)(epsilon=-1.0, lam=model_lam, kappa=1.0, measure_gauss_law=measure_gauss_law) + else: + model = getattr(models, model_name)(alpha=1.0, epsilon=-1.0) + model.kinetic_ions.set_markers( + loading_params=LoadingParameters(ppc=40, seed=1234), + weights_params=WeightsParameters(control_variate=control_variate), + ) + model.kinetic_ions.var.add_background(maxwellians.Maxwellian3D(n=(1.0, None))) + model.kinetic_ions.var.add_initial_condition( + maxwellians.Maxwellian3D(n=(1.0, perturbations.ModesCos(ls=(1,), amps=(1e-1,)))) + ) + model.em_fields.b_field.add_perturbation( + perturbations.ModesSin(ls=(1,), ms=(1,), given_in_basis="2", comp=2, amps=(1e-1,)) + ) + + solver_params = SolverParameters(tol=1e-13, maxiter=3000) + if model_name == "VlasovMaxwellQuasiNeutral": + prop = model.propagators.coupling + prop.options = prop.Options(solver_params=solver_params, **options) + else: + for prop in (model.propagators.maxwell, model.propagators.coupling_va): + prop.options = prop.Options(solver_params=solver_params) + + Simulation( + model=model, + domain=domains.Cuboid(r1=4.0, r2=2.0), + grid=grids.TensorProductGrid(num_elements=(8, 4, 1)), + derham_opts=DerhamOptions(degree=(2, 2, 1)), + ).allocate() + return model + + +def _state(model): + e = model.em_fields.e_field.spline.vector.copy() + b = model.em_fields.b_field.spline.vector.copy() + v = model.kinetic_ions.var.particles.velocities.copy() + return e, b, v + + +def _norm(x): + return x.inner(x) ** 0.5 + + +def _vnorm(v): + """Global 2-norm of the local marker velocities (a collective operation).""" + from maybempi import MPI + + from struphy.propagators.base import Propagator + + val = xp.array([float(xp.sum(v**2))]) + comm = Propagator.derham.comm + if comm is not None and comm.Get_size() > 1: + comm.Allreduce(MPI.IN_PLACE, val, op=MPI.SUM) + return float(val[0]) ** 0.5 + + +# explicit coefficients with c1 != |c2| != |c3| (the species ones are c1 = 1, c2 = c3 = -1 here) +COEFFS = {"species": {}, "explicit": {"c1": 2.0, "c2": -0.5, "c3": -3.0}} + + +@pytest.mark.parametrize("coeffs", list(COEFFS)) +@pytest.mark.parametrize("lam", [1.0, 0.0]) +@pytest.mark.parametrize("precond", ["StiffnessPreconditioner", "MassMatrixPreconditioner"]) +def test_energy_conservation(lam, precond, coeffs): + """The coupling step conserves its energy H exactly (up to the solver tolerance), also for dt >> 1/omega_p and + in the quasi-neutral limit lam=0 (without the initial projection, after which the lam=0 step at frozen + positions exchanges no energy).""" + model = _model(lam=lam, precond=precond, initial_projection=False, **COEFFS[coeffs]) + prop = model.propagators.coupling + b = model.em_fields.b_field.spline.vector + energy_0 = prop.energy() + energy_b_0 = 0.5 * b.inner(prop._M2.dot(b)) + for _ in range(3): + prop(2.0) + assert prop._solver._info["success"] + # H is dominated by the thermal energy: compare its change also with the energy exchanged with the B field + exchanged = abs(0.5 * b.inner(prop._M2.dot(b)) - energy_b_0) + assert exchanged > energy_b_0 + assert abs(prop.energy() - energy_0) < 1e-9 * exchanged + + +@pytest.mark.parametrize("matrix_free", [False, True]) +@pytest.mark.parametrize("lam", [1.0, 0.0]) +def test_discrete_equations(lam, matrix_free): + """The new state satisfies the Crank-Nicolson equations: Ampère's law + lam^2 M1 (e^{n+1} - e^n) / dt = c1 C^T M2 b_mid - c2 j(V_mid) (with V^{n+1} from the pusher, j linear in V at + fixed positions) and Faraday's law (b^{n+1} - b^n) / dt = -C e_mid, with e_mid = e^{n+1} for lam = 0 (without + the initial projection, so that the step is not trivial).""" + dt = 0.5 + model = _model(lam=lam, matrix_free=matrix_free, initial_projection=False, **COEFFS["explicit"]) + prop = model.propagators.coupling + e0, b0, _ = _state(model) + prop._accum(*prop._args_accum) + j0 = prop._accum.vectors[0].copy() + prop(dt) + e1, b1, _ = _state(model) + prop._accum(*prop._args_accum) + j_mid = (j0 + prop._accum.vectors[0]) * 0.5 + b_mid = (b0 + b1) * 0.5 + e_mid = (e0 + e1) * 0.5 if lam > 0.0 else e1 + + curl_b = prop._curl.T.dot(prop._M2.dot(b_mid)) * prop._c1 + current = j_mid * prop._c2 + residual = curl_b - current - prop._M1.dot(e1 - e0) * (lam**2 / dt) + assert _norm(residual) < 1e-8 * (_norm(curl_b) + _norm(current)) + + faraday = (b1 - b0) * (1 / dt) + prop._curl.dot(e_mid) + assert _norm(faraday) < 1e-10 * _norm(b1 - b0) / dt + + if lam == 0.0: + # weakly divergence-free mid-point current + div_j = prop.derham.grad.T.dot(j_mid) + assert _norm(div_j) < 1e-8 * _norm(j_mid) + + +def test_dt_change(): + """The Schur solver and its preconditioner are rebuilt only when dt changes (here with a given n_ref, so the + Schur operator is built before any accumulation); the energy stays conserved across the change.""" + model = _model(n_ref=1.0) + prop = model.propagators.coupling + energy_0 = prop.energy() + prop(0.5) + solver = prop._solver + prop(0.5) + assert prop._solver is solver + prop(1.0) + assert prop._solver is not solver + assert prop._schur_dt == 1.0 + assert abs(prop.energy() - energy_0) < 1e-12 * energy_0 + + +def test_agrees_with_split_propagators(): + """For lam=1 and the species coefficients, one step agrees with MaxwellWeakAmpere followed by + VlasovAmpereCoupling up to the splitting error O(dt^2).""" + diffs = [] + for dt in (0.1, 0.05): + coupled = _model() + coupled.propagators.coupling(dt) + split = _model("VlasovMaxwellOneSpecies") + split.propagators.maxwell(dt) + split.propagators.coupling_va(dt) + # the two simulations have different Derham objects: compare the arrays (global 2-norms) + diff = 0.0 + for x_c, x_s in zip(_state(coupled), _state(split)): + x_c = x_c if isinstance(x_c, xp.ndarray) else x_c.toarray() + x_s = x_s if isinstance(x_s, xp.ndarray) else x_s.toarray() + diff += _vnorm(x_c - x_s) / _vnorm(x_s) + diffs.append(diff) + # local splitting error O(dt^2): halving dt divides the difference by about 4 + assert 3.0 < diffs[0] / diffs[1] < 5.0, diffs + + +def test_quasi_neutral_limit(): + """b and V converge like lam^2 to the lam=0 solution (of the step without the initial projection); the + iteration counts stay bounded as lam -> 0.""" + dt = 0.5 + results = {} + niter = {} + for lam in (1e-1, 1e-2, 0.0): + model = _model(lam=lam, initial_projection=False) + # the same initial data for all lam: the initial Poisson field (Gauss' law lam^2 div E = c2 rho) scales + # like 1 / lam^2, and lam^2 M1 e^n would not vanish + model.em_fields.e_field.spline.vector *= 0.0 + prop = model.propagators.coupling + prop(dt) + niter[lam] = prop._solver._info["niter"] + _, b1, v1 = _state(model) + results[lam] = (b1.toarray(), v1) + + b_ref, v_ref = results[0.0] + errs = [_vnorm(results[k][0] - b_ref) + _vnorm(results[k][1] - v_ref) for k in (1e-1, 1e-2)] + assert errs[1] < 3e-2 * errs[0], errs + assert max(niter.values()) <= 2 * min(niter.values()), niter + + +def _constraint_residual(prop): + """r = c1 C^T M2 b - c2 j(V) at the current markers.""" + b = prop.variables.b.spline.vector + prop._accum(*prop._args_accum) + return prop._curl.T.dot(prop._M2.dot(b)) * prop._c1 - prop._accum.vectors[0] * prop._c2 + + +@pytest.mark.parametrize("matrix_free", [False, True]) +def test_initial_projection(matrix_free): + """lam=0: without the initial projection, the constraint residual r = c1 C^T M2 b - c2 j(V) changes sign in + every step (period-2 oscillation). The projection in the setup makes r vanish, does not depend on dt, keeps e, + and the steps then leave b and V unchanged (frozen positions).""" + # without projection: r^{n+1} = -r^n + model = _model(lam=0.0, initial_projection=False, matrix_free=matrix_free) + prop = model.propagators.coupling + r0 = _constraint_residual(prop) + prop(0.5) + r1 = _constraint_residual(prop) + assert _norm(r0) > 0.0 + assert _norm(r1 + r0) < 1e-10 * _norm(r0) + + # with projection (default), for two dt + projected = [] + for dt in (0.5, 2.0): + model = _model(lam=0.0, matrix_free=matrix_free) + prop = model.propagators.coupling + # accumulate at the initial markers for the scale of the current + prop._accum(*prop._args_accum) + j_scale = abs(prop._c2) * _norm(prop._accum.vectors[0]) + e_init = model.em_fields.e_field.spline.vector.copy() + prop.prepare_time_step(dt) + assert _norm(_constraint_residual(prop)) < 1e-10 * j_scale + e, b, v = _state(model) + # the electric field is not changed by the projection + assert _norm(e - e_init) == 0.0 + projected.append((b.toarray(), v)) + for _ in range(2): + prop(dt) + _, b2, v2 = _state(model) + assert _norm(b2 - b) < 1e-10 * _norm(b) + assert _vnorm(v2 - v) < 1e-10 * _vnorm(v) + assert _norm(_constraint_residual(prop)) < 1e-10 * j_scale + (b_a, v_a), (b_b, v_b) = projected + assert _vnorm(b_a - b_b) < 1e-10 * _vnorm(b_a) + assert _vnorm(v_a - v_b) < 1e-10 * _vnorm(v_a) + + +def test_control_variate_quasi_neutral_rejected(): + """The quasi-neutral limit needs full-f weights.""" + with pytest.raises(AssertionError, match="quasi-neutral"): + _model(control_variate=True, lam=0.0) + + +@pytest.mark.mpi_skip +def test_options(): + """Invalid coefficients are rejected; solver="cg" takes no preconditioner.""" + from struphy.propagators.vlasov_maxwell_coupling import VlasovMaxwellCoupling + + Options = VlasovMaxwellCoupling.Options + with pytest.raises(AssertionError): + Options(lam=-1.0) + with pytest.raises(AssertionError): + Options(c1=0.0) + with pytest.raises(AssertionError): + Options(c2=1.0, c3=-1.0) + with pytest.raises(AssertionError): + Options(precond="HiptmairXu") + with pytest.raises(AssertionError): + Options(matrix_free=1) + assert Options().matrix_free is False + with pytest.raises(AssertionError): + Options(initial_projection=None) + assert Options().initial_projection is True + assert Options().precond == "StiffnessPreconditioner" + assert Options(solver="cg").precond is None + + +@pytest.mark.mpi_skip +def test_equation_parameters(): + """lam = Omega_c x / c and kappa = Omega_p x / c of the species: kappa = alpha lam; for the cyclotron velocity scale + (VlasovMaxwellQuasiNeutral) epsilon = 1, for the light velocity scale (VlasovMaxwellOneSpecies) epsilon lam = 1 and + kappa^2 epsilon = alpha^2 / epsilon, so that the propagator coefficients agree for both normalizations.""" + import math + + from struphy import BaseUnits + from struphy.models import VlasovMaxwellOneSpecies, VlasovMaxwellQuasiNeutral + + base_units = BaseUnits(x=0.01, B=0.5, n=2.0) + for model in (VlasovMaxwellQuasiNeutral(base_units=base_units), VlasovMaxwellOneSpecies(base_units=base_units)): + params = model.kinetic_ions.equation_params + assert math.isclose(params.kappa, params.alpha * params.lam, rel_tol=1e-12) + if isinstance(model, VlasovMaxwellQuasiNeutral): + assert math.isclose(params.epsilon, 1.0, rel_tol=1e-12) + else: + assert math.isclose(params.epsilon * params.lam, 1.0, rel_tol=1e-12) + assert math.isclose(params.kappa**2 * params.epsilon, params.alpha**2 / params.epsilon, rel_tol=1e-12) + + # overrides are passed to the species + params = VlasovMaxwellQuasiNeutral(lam=0.0, kappa=2.0).kinetic_ions.equation_params + assert params.lam == 0.0 and params.kappa == 2.0 + + +def test_model_coefficients(): + """The default propagator coefficients come from the species: lam^2 = (epsilon lam_s)^2, c2 = kappa_s^2 epsilon, + c3 = 1 / epsilon; lam=0 in the model is the quasi-neutral limit (no initial Poisson solve, gauss_error NaN).""" + import math + + model = _model(model_lam=0.5, lam=None, measure_gauss_law=True) + prop = model.propagators.coupling + epsilon = model.kinetic_ions.equation_params.epsilon + assert math.isclose(prop._lam2, (epsilon * 0.5) ** 2) + assert math.isclose(prop._c2, 1.0**2 * epsilon) + assert math.isclose(prop._c3, 1 / epsilon) + # Gauss' law of the initial Poisson solve: lam^2 div E = c2 rho + assert math.isclose(model.initial_poisson.coeffs[0], prop._c2 / prop._lam2) + assert not math.isnan(model.calculate_gauss_error()) + + model = _model(model_lam=0.0, lam=None, measure_gauss_law=True) + assert model.propagators.coupling._lam2 == 0.0 + assert math.isnan(model.calculate_gauss_error()) + assert _norm(model.em_fields.e_field.spline.vector) == 0.0 diff --git a/src/struphy/propagators/vlasov_maxwell_coupling.py b/src/struphy/propagators/vlasov_maxwell_coupling.py new file mode 100644 index 000000000..18234cd8c --- /dev/null +++ b/src/struphy/propagators/vlasov_maxwell_coupling.py @@ -0,0 +1,569 @@ +"Particle and FEEC variables are updated." + +import logging +from dataclasses import dataclass + +import cunumpy as xp +from feectools.linalg.solvers import inverse +from line_profiler import profile +from maybempi import MPI + +from struphy.feec import preconditioner +from struphy.io.options import LiteralOptions, OptionsBase +from struphy.linear_algebra.solver import SolverParameters +from struphy.models.variables import FEECVariable, PICVariable +from struphy.pic.accumulation.kernels.marker_operator_symm_dot import marker_operator_symm_dot +from struphy.pic.accumulation.kernels.vlasov_maxwell import vlasov_maxwell +from struphy.pic.accumulation.kernels.vlasov_maxwell_mf import vlasov_maxwell_mf +from struphy.pic.accumulation.particles_to_grid import Accumulator, AccumulatorVector, MarkerOperator +from struphy.pic.pushing.kernels.push_v_with_efield import push_v_with_efield +from struphy.pic.pushing.pusher import Pusher +from struphy.propagators.base import Propagator +from struphy.utils.utils import check_option + +logger = logging.getLogger("struphy") + + +class VlasovMaxwellCoupling(Propagator): + r"""PIC-FEEC discretization of the following equations: + find :math:`\mathbf E \in H(\textnormal{curl})`, :math:`\mathbf B \in H(\textnormal{div})` and :math:`f` such that + + .. math:: + + \lambda^2 \int_{\Omega} \frac{\partial \mathbf E}{\partial t} \cdot \mathbf F \,\textrm d \mathbf x + &= c_1 \int_{\Omega} \mathbf B \cdot \nabla \times \mathbf F \,\textrm d \mathbf x + - c_2 \int_{\Omega} \int_{\mathbb{R}^3} f \mathbf{v} \cdot \mathbf{F} \, \text{d}^3 \mathbf{v} \,\textrm d \mathbf x + \qquad \forall \, \mathbf F \in H(\textnormal{curl})\,, + \\[2mm] + \frac{\partial \mathbf B}{\partial t} + \nabla\times\mathbf E &= 0\,, + \\[2mm] + \frac{\partial f}{\partial t} + c_3\, \mathbf{E} \cdot \frac{\partial f}{\partial \mathbf{v}} &= 0\,, + + with :math:`\lambda \geq 0`, :math:`c_1 > 0` and :math:`c_2 c_3 > 0`. It solves the steps of + :class:`~struphy.propagators.maxwell_weak_ampere.MaxwellWeakAmpere` and + :class:`~struphy.propagators.vlasov_ampere_coupling.VlasovAmpereCoupling` in one implicit step (for + :math:`\lambda = 1`, :math:`c_1 = 1`, :math:`c_2 = \alpha^2/\varepsilon`, :math:`c_3 = 1/\varepsilon`, the light + velocity scale). The defaults are the coefficients of the species' :ref:`equation parameters `, + :math:`\lambda = \varepsilon \lambda_\textrm{s}`, :math:`c_2 = \kappa_\textrm{s}^2 \varepsilon`, + :math:`c_3 = 1/\varepsilon`, i.e. :math:`\lambda = \lambda_\textrm{s}`, :math:`c_2 = \kappa_\textrm{s}^2`, + :math:`c_3 = 1` for the cyclotron velocity scale. + For :math:`\lambda = 0` the displacement current vanishes, :math:`c_1 \nabla \times \mathbf B = c_2 \mathbf J` + (quasi-neutral limit). + + With the FE coefficients :math:`\mathbf e \in \mathbb R^{N_1}` of :math:`\mathbf E \in V^1_h` and + :math:`\mathbf b \in \mathbb R^{N_2}` of :math:`\mathbf B \in V^2_h`, and :math:`f` approximated by :math:`N` + markers with positions :math:`\boldsymbol \eta_p` (logical, fixed in this step), velocities + :math:`\mathbf V = (\mathbf v_p)_p \in \mathbb R^{3N}` and weights :math:`\mathbb W = \textnormal{diag}(w_p)`, + the semi-discrete system reads + + .. math:: + + \lambda^2\, \mathbb M_1 \dot{\mathbf e} = c_1\, \mathbb C^\top \mathbb M_2 \mathbf b + - c_2\, (\mathbf \Lambda^1)^\top \, \overline{DF^{-1}} \, \mathbb W \, \mathbf V \,, + \qquad + \dot{\mathbf b} = - \mathbb C \mathbf e\,, + \qquad + \dot{\mathbf V} = c_3\, \overline{DF^{-\top}} \, \mathbf \Lambda^1 \, \mathbf e\,, + + with the notation of :class:`~struphy.propagators.vlasov_ampere_coupling.VlasovAmpereCoupling`. It conserves the + energy (:meth:`~struphy.propagators.vlasov_maxwell_coupling.VlasovMaxwellCoupling.energy`) + + .. math:: + + H = \frac{\lambda^2}{2}\, \mathbf e^\top \mathbb M_1 \mathbf e + \frac{c_1}{2}\, \mathbf b^\top \mathbb M_2 \mathbf b + + \frac{c_2}{2 c_3}\, \mathbf V^\top \mathbb W \mathbf V\,. + + :ref:`time_discret`: Crank-Nicolson (implicit mid-point), with + :math:`\mathbf e^\textnormal{mid} = \frac 12 (\mathbf e^{n+1} + \mathbf e^n)` etc., + + .. math:: + + \lambda^2\, \mathbb M_1 \frac{\mathbf e^{n+1} - \mathbf e^n}{\Delta t} + = c_1\, \mathbb C^\top \mathbb M_2 \mathbf b^\textnormal{mid} + - c_2\, (\mathbf \Lambda^1)^\top \, \overline{DF^{-1}} \, \mathbb W \, \mathbf V^\textnormal{mid}\,, + \qquad + \frac{\mathbf b^{n+1} - \mathbf b^n}{\Delta t} = - \mathbb C \mathbf e^\textnormal{mid}\,, + \qquad + \frac{\mathbf V^{n+1} - \mathbf V^n}{\Delta t} = c_3\, \overline{DF^{-\top}} \, \mathbf \Lambda^1 \, \mathbf e^\textnormal{mid}\,, + + which conserves :math:`H` exactly (up to the solver tolerance; for constant weights, i.e. without control + variate). Eliminating :math:`\mathbf b^\textnormal{mid}` and :math:`\mathbf V^\textnormal{mid}` gives the Schur + system for the mid-point field + + .. math:: + + \left(\lambda^2\, \mathbb M_1 + \frac{\Delta t^2}{4}\, \mathbb K\right) \mathbf e^\textnormal{mid} + = \lambda^2\, \mathbb M_1 \mathbf e^n + + \frac{\Delta t}{2} \left( c_1\, \mathbb C^\top \mathbb M_2 \mathbf b^n - c_2\, \mathbf j \right)\,, + \qquad + \mathbb K = c_1\, \mathbb C^\top \mathbb M_2 \mathbb C + c_2 c_3\, \mathbb A\,, + + with the accumulated matrix :math:`\mathbb A` and vector :math:`\mathbf j` of + :class:`~struphy.propagators.vlasov_ampere_coupling.VlasovAmpereCoupling` (kernel + :func:`~struphy.pic.accumulation.kernels.vlasov_maxwell.vlasov_maxwell_kernels.vlasov_maxwell`). It is solved + with ``solver`` and ``precond``. Then + + .. math:: + + \mathbf b^{n+1} = \mathbf b^n - \Delta t\, \mathbb C \mathbf e^\textnormal{mid}\,, + \qquad + \mathbf e^{n+1} = 2 \mathbf e^\textnormal{mid} - \mathbf e^n \quad (\lambda > 0)\,, + + and the velocities are updated by + :func:`~struphy.pic.pushing.kernels.push_v_with_efield.push_v_with_efield_kernels.push_v_with_efield` with + :math:`\mathbf e^\textnormal{mid}`. + + **Quasi-neutral limit** :math:`\lambda = 0`. The curl-curl operator :math:`\mathbb C^\top \mathbb M_2 \mathbb C` + is singular on the discrete gradients (and harmonic fields); this kernel is absorbed by the plasma term + :math:`c_2 c_3 \mathbb A`, a particle approximation of the density-weighted mass matrix :math:`\mathbb M_1`. + The Schur operator is therefore positive definite for positive weights (full-f, no control variate) and plasma + everywhere in the domain. The scheme then imposes + :math:`c_1\, \mathbb C^\top \mathbb M_2 \mathbf b^\textnormal{mid} = c_2\, (\mathbf \Lambda^1)^\top + \overline{DF^{-1}} \mathbb W \mathbf V^\textnormal{mid}`, hence a weakly divergence-free mid-point current + (multiply by :math:`\mathbb G^\top`); the gradient part of :math:`\mathbf E` is its Lagrange multiplier. + :math:`\mathbf e^n` drops out of the system and only :math:`\mathbf e^\textnormal{mid}` is determined: the + propagator stores :math:`\mathbf e^{n+1} = \mathbf e^\textnormal{mid}` (the Crank-Nicolson update + :math:`2 \mathbf e^\textnormal{mid} - \mathbf e^n` would carry any error of :math:`\mathbf e^n` as an undamped + even-odd oscillation). The energy :math:`H` (without the electric part) is still conserved exactly. + + The constraint holds at the mid-point only: its residual + :math:`\mathbf r^n = c_1 \mathbb C^\top \mathbb M_2 \mathbf b^n - c_2 \mathbf j(\mathbf V^n)` (linear in + :math:`\mathbf V` at fixed positions) satisfies :math:`\mathbf r^{n+1} = - \mathbf r^n`. Initial data that violate + the constraint (e.g. by marker noise in :math:`\mathbf J`) would thus oscillate with period 2 without damping. With + ``initial_projection`` (default), the setup replaces :math:`(\mathbf b^0, \mathbf V^0)` by the mid-point state of + one step, :math:`\mathbf b^0 - \frac{\Delta t}{2} \mathbb C \mathbf e^\textnormal{mid}` and + :math:`\mathbf V^0 + \frac{\Delta t}{2} c_3 \overline{DF^{-\top}} \mathbf \Lambda^1 \mathbf e^\textnormal{mid}` + (an implicit Euler step of :math:`\Delta t / 2`), which satisfies the constraint exactly and does not depend on + :math:`\Delta t`: a projection onto the constraint. The particle pushers of the splitting change + :math:`\mathbf r` again by :math:`O(\Delta t)` in every step. + + **Preconditioner.** With :math:`\mathbb A \approx \bar n\, \mathbb M_1`, the Schur operator is approximately the + stabilized curl-curl operator + :math:`\mathbb C^\top \left(\frac{\Delta t^2 c_1}{4} \mathbb M_2\right) \mathbb C + + \left(\lambda^2 + \frac{\Delta t^2 c_2 c_3 \bar n}{4}\right) \mathbb M_1`, preconditioned (default) by + :class:`~struphy.feec.preconditioner.StiffnessPreconditioner`, built for the current time step + (:meth:`~struphy.propagators.vlasov_maxwell_coupling.VlasovMaxwellCoupling.prepare_time_step`). Its + stabilization stays positive for :math:`\lambda = 0`. + """ + + class Variables: + """Container for variables advanced by :class:`~struphy.propagators.vlasov_maxwell_coupling.VlasovMaxwellCoupling`. + + Attributes + ---------- + e : FEECVariable + Electric-field variable in ``"Hcurl"`` space. + b : FEECVariable + Magnetic-field variable in ``"Hdiv"`` space. + ions : PICVariable + Particle variable in ``"Particles6D"`` space. + """ + + def __init__(self): + self._e: FEECVariable = None + self._b: FEECVariable = None + self._ions: PICVariable = None + + @property + def e(self) -> FEECVariable: + return self._e + + @e.setter + def e(self, new): + assert isinstance(new, FEECVariable) + assert new.space == "Hcurl" + self._e = new + + @property + def b(self) -> FEECVariable: + return self._b + + @b.setter + def b(self, new): + assert isinstance(new, FEECVariable) + assert new.space == "Hdiv" + self._b = new + + @property + def ions(self) -> PICVariable: + return self._ions + + @ions.setter + def ions(self, new): + assert isinstance(new, PICVariable) + assert new.space == "Particles6D" + self._ions = new + + def __init__(self): + self.variables = self.Variables() + + @dataclass(repr=False) + class Options(OptionsBase): + r"""Configuration options for :class:`~struphy.propagators.vlasov_maxwell_coupling.VlasovMaxwellCoupling`. + + Parameters + ---------- + lam : float, default=None + :math:`\lambda \geq 0`, with the coefficient :math:`\lambda^2` of the displacement current + :math:`\partial_t \mathbf E`; ``lam=0.0`` is the quasi-neutral limit (needs full-f weights, no control + variate). If ``None``, :math:`\varepsilon \lambda_\textrm{s}` of the species (``lam`` and ``epsilon`` of + its equation parameters). + + c1 : float, default=1.0 + Coefficient :math:`c_1 > 0` of :math:`\nabla \times \mathbf B` in Ampère's law. + + c2 : float, default=None + Coefficient of the current density in Ampère's law. If ``None``, :math:`\kappa_\textrm{s}^2 \varepsilon` of + the species (``kappa`` and ``epsilon`` of its equation parameters). + + c3 : float, default=None + Coefficient of the electric force in the Vlasov equation (:math:`c_2 c_3 > 0`). If ``None``, + :math:`1/\varepsilon` of the species. + + solver : LiteralOptions.OptsSymmSolver, default="pcg" + Symmetric iterative solver of the Schur system, passed to ``feectools.linalg.solvers.inverse``. + + precond : LiteralOptions.OptsMaxwellPrecond, default="StiffnessPreconditioner" + Preconditioner of the Schur system. + + - ``"StiffnessPreconditioner"``: :class:`~struphy.feec.preconditioner.StiffnessPreconditioner` with + ``derivative="curl"`` for the stabilized curl-curl approximation of the Schur operator (see the class + docstring), rebuilt when the time step changes. + - ``"MassMatrixPreconditioner"``, ``"PolarMassMatrixPreconditioner"``: preconditioner of + :math:`\mathbb M_1`, sufficient for :math:`\Delta t \lesssim h` (polar splines are not supported by + this propagator yet). + - ``None``: no preconditioner (set automatically, with a warning, for ``solver="cg"``). + + solver_params : SolverParameters, default=None + Solver controls. If ``None``, defaults to ``SolverParameters()``. + + n_ref : float, default=None + Reference density :math:`\bar n` in the stiffness preconditioner. If ``None``, the Rayleigh quotient + :math:`\mathbf 1^\top \mathbb A \mathbf 1 / \mathbf 1^\top \mathbb M_1 \mathbf 1` of the accumulated + matrix at the setup (a density average). + + matrix_free : bool, default=False + Apply the particle matrix :math:`\mathbb A` matrix-free, with the filling + :math:`w_p G^{-1}(\boldsymbol \eta_p)` stored per marker + (:class:`~struphy.pic.accumulation.particles_to_grid.MarkerOperator`), instead of assembling it: one spline + evaluation and one deposition per marker and solver iteration instead of the products of all pairs of + basis functions per marker. Needs 6 auxiliary marker columns. Ignored on the CuPy backend. + + initial_projection : bool, default=True + Only for ``lam=0``: in the setup (first :meth:`prepare_time_step`), replace :math:`(\mathbf b, \mathbf V)` + by the mid-point state of one step, which satisfies :math:`c_1 \nabla \times \mathbf B = c_2 \mathbf J` + exactly (see the class docstring). Without it, initial data that violate the constraint lead to an undamped + period-2 oscillation of :math:`\mathbf b` and :math:`\mathbf V`. Not applied to a restored (restart) state. + """ + + lam: float = None + c1: float = 1.0 + c2: float = None + c3: float = None + solver: LiteralOptions.OptsSymmSolver = "pcg" + precond: LiteralOptions.OptsMaxwellPrecond = "StiffnessPreconditioner" + solver_params: SolverParameters = None + n_ref: float = None + matrix_free: bool = False + initial_projection: bool = True + + def __post_init__(self): + # checks + check_option(self.solver, LiteralOptions.OptsSymmSolver) + check_option(self.precond, LiteralOptions.OptsMaxwellPrecond) + if self.lam is not None: + assert self.lam >= 0.0, f"lam must be >= 0, got {self.lam}." + assert self.c1 > 0.0, f"c1 must be > 0, got {self.c1}." + if self.c2 is not None and self.c3 is not None: + assert self.c2 * self.c3 > 0.0, f"c2 * c3 must be > 0, got c2={self.c2}, c3={self.c3}." + if self.n_ref is not None: + assert self.n_ref > 0.0, f"n_ref must be > 0, got {self.n_ref}." + assert isinstance(self.matrix_free, bool) + assert isinstance(self.initial_projection, bool) + if self.solver == "cg" and self.precond is not None: + logger.warning(f'solver="cg" takes no preconditioner; ignoring precond="{self.precond}".') + self.precond = None + + # defaults + if self.solver_params is None: + self.solver_params = SolverParameters() + + @property + def options(self) -> Options: + if not hasattr(self, "_options"): + self._options = self.Options() + return self._options + + @options.setter + def options(self, new): + assert isinstance(new, self.Options) + self._options = new + logger.info(f"\nNew options for propagator '{self.__class__.__name__}':\n{self._options}") + + @profile + def allocate(self): + # coefficients + equation_params = self.variables.ions.species.equation_params + epsilon = equation_params.epsilon + # lam^2 = (epsilon lam_s)^2, c2 = kappa_s^2 epsilon, c3 = 1 / epsilon (see the equation parameters) + self._lam2 = (epsilon * equation_params.lam) ** 2 if self.options.lam is None else self.options.lam**2 + self._c1 = self.options.c1 + self._c2 = equation_params.kappa**2 * epsilon if self.options.c2 is None else self.options.c2 + self._c3 = 1.0 / epsilon if self.options.c3 is None else self.options.c3 + assert self._c2 * self._c3 > 0.0, f"c2 * c3 must be > 0, got c2={self._c2}, c3={self._c3}." + + assert not self.derham.polar_splines, "VlasovMaxwellCoupling does not support polar splines yet." + + particles = self.variables.ions.particles + if particles.control_variate: + assert self._lam2 > 0.0, ( + "The quasi-neutral limit lam=0 needs full-f weights: with the control variate the accumulated " + "matrix is indefinite and the Schur system is not well-posed." + ) + logger.warning( + "VlasovMaxwellCoupling with control variate: the energy is not conserved exactly (the weights change)." + ) + + self._info = self.options.solver_params.info + + # accumulation of the matrix A (assembled, or matrix-free) and the vector j; + # the matrix-free kernels have no CUDA version + self._matrix_free = self.options.matrix_free and xp.get_backend() != "cupy" + if self._matrix_free: + # the vector, and the filling of the matrix per marker (in the first free columns of the marker array) + column = particles.args_markers.first_free_idx + assert particles.n_cols_aux >= 6, "matrix_free needs 6 free (auxiliary) marker columns." + self._args_accum = (column,) + self._accum = AccumulatorVector( + particles, + "Hcurl", + vlasov_maxwell_mf, + self.mass_ops, + self.domain.args_domain, + ) + self._A_particles = MarkerOperator( + particles, + "Hcurl", + marker_operator_symm_dot, + self.mass_ops, + self.domain.args_domain, + (1, column), + symmetry="symm", + ) + else: + self._args_accum = () + self._accum = Accumulator( + particles, + "Hcurl", + vlasov_maxwell, + self.mass_ops, + self.domain.args_domain, + add_vector=True, + symmetry="symm", + ) + self._A_particles = self._accum.operators[0] + + # operators + self._M1 = self.mass_ops.M1 + self._M2 = self.mass_ops.M2 + self._curl = self.derham.curl + self._curl_curl = self._curl.T @ self._M2 @ self._curl + + # reference density for the stiffness preconditioner: density average of the accumulated matrix + if self.options.n_ref is None: + self._accum(*self._args_accum) + ones = self.derham.V1.zeros() + for block in ones.blocks: + block._data[:] = 1.0 + ones_A = ones.inner(self._A_particles.dot(ones)) + ones_M1 = ones.inner(self._M1.dot(ones)) + if ones_A > 0.0: + self._n_ref = ones_A / ones_M1 + else: + self._n_ref = 1.0 + logger.warning("VlasovMaxwellCoupling: no positive density average of the markers, using n_ref=1.") + else: + self._n_ref = self.options.n_ref + + # preconditioner (the stiffness preconditioner depends on dt and is built in prepare_time_step) + self._precond = self.options.precond + if self._precond in (None, "StiffnessPreconditioner"): + self._mass_pc = None + else: + pc_class = getattr(preconditioner, self._precond) + self._mass_pc = pc_class(self._M1) + + # solver kwargs (the Schur operator and its solver are built in prepare_time_step) + kwargs = self.options.solver_params.__dict__.copy() + kwargs.pop("info") + self._solver_kwargs = kwargs + self._solver = None + self._schur_dt = None + self._projected = False + + # buffers (the pusher reads e_mid) + self._e_mid = self.derham.V1.zeros() + self._e_new = self.derham.V1.zeros() + self._e_rhs = self.derham.V1.zeros() + self._e_tmp = self.derham.V1.zeros() + self._b_new = self.derham.V2.zeros() + self._b_tmp = self.derham.V2.zeros() + + # particle pusher: v^{n+1} = v^n + dt * c3 * DF^{-T} Lambda^1 e_mid + args_kernel = ( + self.derham.args_derham, + self._e_mid.blocks[0]._data, + self._e_mid.blocks[1]._data, + self._e_mid.blocks[2]._data, + self._c3, + ) + + self._pusher = Pusher( + particles, + push_v_with_efield, + args_kernel, + self.domain.args_domain, + alpha_in_kernel=1.0, + pushes_eta=False, + ) + + def prepare_time_step(self, dt: float): + """Schur operator :math:`\\lambda^2 \\mathbb M_1 + \\frac{\\Delta t^2}{4} \\mathbb K`, its preconditioner and + solver for the time step ``dt``: built in the setup, rebuilt when dt changes. The particle matrix enters by + reference (assembled, or matrix-free reading the fillings stored in the markers), so it is updated by the + accumulation in every step. For ``lam=0`` (and ``initial_projection``), the first call also projects the + state onto the constraint :math:`c_1 \\mathbb C^\\top \\mathbb M_2 \\mathbf b = c_2 \\mathbf j`.""" + if dt == self._schur_dt: + return + + schur = (dt**2 * self._c1 / 4) * self._curl_curl + (dt**2 * self._c2 * self._c3 / 4) * self._A_particles + if self._lam2 > 0.0: + schur = self._lam2 * self._M1 + schur + + if self._precond == "StiffnessPreconditioner": + pc = preconditioner.StiffnessPreconditioner( + self.mass_ops, + "curl", + sigma=self._lam2 + dt**2 * self._c2 * self._c3 * self._n_ref / 4, + mid=(dt**2 * self._c1 / 4) * self._M2, + stab="M1", + ) + else: + pc = self._mass_pc + + kwargs = dict(self._solver_kwargs) + if pc is not None: + kwargs["pc"] = pc + self._solver = inverse(schur, self.options.solver, **kwargs) + self._schur_dt = dt + + # lam = 0: consistent initial state, once (in the setup, before the initial data are saved) + if self._lam2 == 0.0 and self.options.initial_projection and not self._projected: + self._project_onto_constraint(dt) + self._projected = True + + def energies(self) -> tuple[float, float, float]: + r"""The three parts :math:`\frac{\lambda^2}{2} \mathbf e^\top \mathbb M_1 \mathbf e`, + :math:`\frac{c_1}{2} \mathbf b^\top \mathbb M_2 \mathbf b` and + :math:`\frac{c_2}{2 c_3} \mathbf V^\top \mathbb W \mathbf V` of the energy :meth:`energy` (summed over all + MPI ranks; a collective operation).""" + e = self.variables.e.spline.vector + b = self.variables.b.spline.vector + particles = self.variables.ions.particles + + kinetic = xp.array([0.5 * xp.sum(particles.weights * xp.sum(particles.velocities**2, axis=0))]) + if self.derham.comm is not None: + self.derham.comm.Allreduce(MPI.IN_PLACE, kinetic, op=MPI.SUM) + if particles.clone_config is not None: + particles.clone_config.inter_comm.Allreduce(MPI.IN_PLACE, kinetic, op=MPI.SUM) + + energy_e = self._lam2 / 2 * e.inner(self._M1.dot(e, out=self._e_tmp)) if self._lam2 > 0.0 else 0.0 + energy_b = self._c1 / 2 * b.inner(self._M2.dot(b, out=self._b_tmp)) + energy_f = self._c2 / self._c3 * float(kinetic[0]) + return energy_e, energy_b, energy_f + + def energy(self) -> float: + r"""The energy :math:`H = \frac{\lambda^2}{2} \mathbf e^\top \mathbb M_1 \mathbf e + \frac{c_1}{2} + \mathbf b^\top \mathbb M_2 \mathbf b + \frac{c_2}{2 c_3} \mathbf V^\top \mathbb W \mathbf V` conserved by the + propagator (summed over all MPI ranks; a collective operation).""" + return sum(self.energies()) + + def _mid_point_field(self, dt: float): + """Accumulate at the current markers and solve the Schur system for the mid-point field (into ``_e_mid``, + with updated ghost regions, which the pusher kernel reads).""" + en = self.variables.e.spline.vector + bn = self.variables.b.spline.vector + + # accumulate A (or its fillings per marker) and j at the current marker positions and velocities + self._accum(*self._args_accum) + + # rhs = lam^2 M1 e^n + dt/2 (c1 C^T M2 b^n - c2 j) + rhs = self._curl.T.dot(self._M2.dot(bn, out=self._b_tmp), out=self._e_rhs) + rhs *= self._c1 + j = self._accum.vectors[0] + j.copy(out=self._e_tmp) + self._e_tmp *= self._c2 + rhs -= self._e_tmp + rhs *= dt / 2 + if self._lam2 > 0.0: + self._M1.dot(en, out=self._e_tmp) + self._e_tmp *= self._lam2 + rhs += self._e_tmp + + e_mid = self._solver.dot(rhs, out=self._e_mid) + e_mid.update_ghost_regions() + return e_mid, self._solver._info + + def _project_onto_constraint(self, dt: float): + """For lam = 0: replace (b, V) by the mid-point state of one step, which satisfies + c1 C^T M2 b = c2 j(V) exactly (independent of dt). The electric field is kept (e.g. from the initial Poisson + solve): e_mid is an impulse of size 1/dt here, and e does not enter the steps for lam = 0.""" + bn = self.variables.b.spline.vector + e_mid, _ = self._mid_point_field(dt) + + # b_mid = b^n - dt/2 C e_mid, V_mid = V^n + dt/2 c3 DF^{-T} Lambda^1 e_mid + b_new = self._curl.dot(e_mid, out=self._b_new) + b_new *= -dt / 2 + b_new += bn + self._pusher(dt / 2) + self.update_feec_variables(b=b_new) + logger.info("VlasovMaxwellCoupling: projected the initial state onto c1 curl B = c2 J (lam = 0).") + + @profile + def __call__(self, dt): + en = self.variables.e.spline.vector + bn = self.variables.b.spline.vector + + # Schur operator for the current time step (no-op unless dt changed since the setup) + self.prepare_time_step(dt) + + # mid-point field + e_mid, info = self._mid_point_field(dt) + + # new b coeffs: b^{n+1} = b^n - dt C e_mid + b_new = self._curl.dot(e_mid, out=self._b_new) + b_new *= -dt + b_new += bn + + # new e coeffs: Crank-Nicolson for lam > 0, the mid-point field for lam = 0 + e_new = e_mid.copy(out=self._e_new) + if self._lam2 > 0.0: + e_new *= 2.0 + e_new -= en + + # update velocities with e_mid + self._pusher(dt) + + # update weights + if self.variables.ions.species.weights_params.control_variate: + self.variables.ions.particles.update_weights() + + # write new coeffs into self.variables + diffs = self.update_feec_variables(e=e_new, b=b_new) + + if self._info and MPI.COMM_WORLD.Get_rank() == 0: + logger.info(f"Status for VlasovMaxwellCoupling: {info['success']}") + logger.info(f"Iterations for VlasovMaxwellCoupling: {info['niter']}") + logger.info(f"Maxdiff e for VlasovMaxwellCoupling: {diffs['e']}") + logger.info(f"Maxdiff b for VlasovMaxwellCoupling: {diffs['b']}") + logger.info("")