diff --git a/doc/sections/subsections/models-all.rst b/doc/sections/subsections/models-all.rst index b896f7082..454465b9c 100644 --- a/doc/sections/subsections/models-all.rst +++ b/doc/sections/subsections/models-all.rst @@ -49,6 +49,11 @@ Gyrokinetic :undoc-members: :show-inheritance: +.. autoclass:: struphy.models.DriftKineticMaxwell + :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.ToyDrift :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/propagators-avail.rst b/doc/sections/subsections/propagators-avail.rst index ac971567a..37c30958c 100644 --- a/doc/sections/subsections/propagators-avail.rst +++ b/doc/sections/subsections/propagators-avail.rst @@ -202,6 +202,11 @@ Vlasov-Maxwell coupling :exclude-members: options, allocate :show-inheritance: +.. automodule:: struphy.propagators.drift_kinetic_maxwell_coupling + :members: + :exclude-members: options, allocate + :show-inheritance: + MHD-kinetic coupling ^^^^^^^^^^^^^^^^^^^^ diff --git a/examples/DriftKineticMaxwell/inertial_alfven_wave/params_inertial_alfven_wave.py b/examples/DriftKineticMaxwell/inertial_alfven_wave/params_inertial_alfven_wave.py new file mode 100644 index 000000000..7d67ebc8a --- /dev/null +++ b/examples/DriftKineticMaxwell/inertial_alfven_wave/params_inertial_alfven_wave.py @@ -0,0 +1,188 @@ +"""Inertial Alfvén wave with the DriftKineticMaxwell model, with and without displacement current (`--lam 1` or +`--lam 0`) and with or without the E x B current (`--exb-current`).""" + +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[_exb]).") +parser.add_argument( + "--lam", + type=float, + default=0.0, + help="lam = v_th / c, coefficient lam^2 of the displacement current (0: quasi-neutral).", +) +parser.add_argument( + "--exb-current", action="store_true", help="Include the E x B current of the guiding centers in Ampere's law." +) +parser.add_argument("--dt", type=float, default=0.5, help="Time step.") +parser.add_argument("--tend", type=float, default=80.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}" + ("_exb" if args.exb_current else "") +if args.test_run: + args.tend = 2 * args.dt + +name = "Inertial Alfven wave in the drift-kinetic Maxwell model" + +description = f""" +Shear (inertial) Alfven wave in a cold guiding-center plasma with the DriftKineticMaxwell model: drift-kinetic +particles (PIC, Particles5D) coupled to Maxwell's equations through the implicit, energy conserving +DriftKineticMaxwellCoupling step, in a periodic slab with the uniform field B0 = e_z. Here lam = {args.lam:g} +(coefficient lam^2 of the displacement current), E x B current {"on" if args.exb_current else "off"}. + +Physics. The model is normalized with the thermal velocity v_th of the species (t = x / v_th, E in units of +v_th B), with the species parameters epsilon = 1 / (Omega_c t), lam_s = Omega_c x / c and kappa_s = Omega_p x / c. +Ampere's law reads + lam^2 dE/dt + chi n0 / |B0|^2 dE_perp/dt = curl B - c2 J_gc, +with the guiding-center current J_gc, lam = epsilon lam_s = v_th / c (here an option, 0 or 1), the polarization +coefficient chi = (epsilon kappa_s)^2 = v_th^2 / v_A^2 and c2 = epsilon kappa_s^2; the parallel force has +c3 = 1 / epsilon. The drift-kinetic ordering removes the polarization drift from the guiding-center motion; the chi +term adds it back. Without it, the quasi-neutral limit lam = 0 would be singular for flute modes. For a cold plasma +and k = (kx, 0, kz), the shear mode (E_x, E_z, B_y) obeys, without the E x B current and for lam = 0, + omega^2 = (1 / chi) kz^2 / (1 + kx^2 d^2), d^2 = 1 / P, P = c2 c3 n0 = kappa_s^2, +the inertial Alfven wave with v_A^2 = 1 / chi and the inertial length d = 1 / kappa_s. For lam > 0, and with the +E x B current (which couples the shear and the fast wave through the Hall term of the single species), the +frequencies are the roots of the determinant of the 3x3 cold-plasma dispersion matrix (computed in the pproc +script). + +Setup. epsilon = 0.5, lam_s = 2, kappa_s = 6 (chi = 9, v_A = 1/3, d = 1/6, c2 = 18, c3 = 2; lam = 1 would mean +v_th = c and only tests the displacement-current branch); box 1 x 1 x 4 with 8 x 1 x 8 cells and cubic splines in +x and z, kx = 2 pi, kz = pi / 2 (kx^2 d^2 = 1.1); cold guiding centers (vth = 1e-6, full-f weights, 50 markers per +cell, seed 1234), initial B_y = 1e-2 sin(kx x + kz z) (2-form component). Time step dt = {args.dt:g} +(omega dt = 0.18) up to t = {args.tend:g}, Lie-Trotter splitting of the coupling step with the parallel and the +E* x b0 guiding-center pushers. For lam = 0 the coupling step first projects the initial state onto the constraint +of the quasi-neutral limit (option initial_projection, default). + +Expected result. Analytical frequencies of the shear mode: omega = 0.3616 (lam = 0), 0.2643 (lam = 0 with the +E x B current) and 0.3428 (lam = 1; the B_y perturbation then also excites the upper branch, omega = 8.70, with a +small amplitude). The pproc fits omega to the magnetic energy en_B ~ cos^2(omega t) and compares. Measured with +this example (serial, 2026-10-23): -1.0 % for lam = 0; with the verification test of the same setup -0.25 % with +the E x B current and -0.7 % for lam = 1 (dt = 0.1). The deviation decreases with the number of markers (-0.35 % +with 200 markers per cell). The total energy en_tot is conserved by the whole splitting up to the solver tolerance +(cold plasma, uniform B0): relative change 5e-6 with the default tolerance 1e-8, 1e-9 with 1e-12. + +References: +- A. J. Brizard and T. S. Hahm, "Foundations of nonlinear gyrokinetic theory", Rev. Mod. Phys. 79, 421 (2007). +- C. K. Birdsall and A. B. Langdon, "Plasma Physics via Computer Simulation" (IOP, 1991). +""" + +# ------------------ +# Import Struphy API +# ------------------ + +import logging + +from struphy import ( + BaseUnits, + DerhamOptions, + EnvironmentOptions, + LoadingParameters, + SavingParameters, + Simulation, + Time, + WeightsParameters, + domains, + equils, + grids, + maxwellians, + perturbations, + set_logging_level, +) + +set_logging_level(logging.WARNING) + +# --------------------- +# Instance of the model +# --------------------- + +from struphy.models import DriftKineticMaxwell + +EPSILON, LAM, KAPPA = 0.5, 2.0, 6.0 +LX, LZ = 1.0, 4.0 +VTH = 1e-6 + +model = DriftKineticMaxwell(base_units=BaseUnits(kBT=1.0), epsilon=EPSILON, lam=LAM, kappa=KAPPA) + +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") +domain = domains.Cuboid(r1=LX, r2=1.0, r3=LZ) +equil = equils.HomogenSlab(B0z=1.0) +grid = grids.TensorProductGrid(num_elements=(8, 1, 8)) +derham_opts = DerhamOptions(degree=(3, 1, 3)) + +sim = Simulation( + model=model, + name=name, + description=description, + params_path=__file__, + env=env, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, +) + +# ------------------- +# Particle parameters +# ------------------- + +model.kinetic_ions.set_markers( + loading_params=LoadingParameters( + ppc=5 if args.test_run else 50, + moments=(0.0, 0.0, VTH, VTH), + B0=1.0, + seed=1234, + ), + weights_params=WeightsParameters(control_variate=False), + saving_params=SavingParameters(n_markers=0), +) + +# ------------------ +# Propagator options +# ------------------ + +prop = model.propagators.coupling +if args.exb_current: + prop.options = prop.Options(lam=args.lam) +else: + # without the E x B current the Schur system is symmetric + prop.options = prop.Options(lam=args.lam, exb_current=False, solver="pcg") +model.propagators.push_para.options = model.propagators.push_para.Options() +# with the electric field (not a potential) only the explicit E* x b0 push is consistent +model.propagators.push_bxe.options = model.propagators.push_bxe.Options(algo="explicit") + +# ------------------ +# Initial conditions +# ------------------ + +cold = dict(n=(1.0, None), vth_para=(VTH, None), vth_perp=(VTH, None), B0=1.0) +model.kinetic_ions.var.add_background(maxwellians.GyroMaxwellian2D(**cold)) +model.em_fields.b_field.add_perturbation( + perturbations.ModesSin(ls=(1,), ns=(1,), given_in_basis="2", comp=1, amps=(1e-2,)) +) + +if __name__ == "__main__": + # post-processing and plots: pproc_inertial_alfven_wave.py + sim.run() diff --git a/examples/DriftKineticMaxwell/inertial_alfven_wave/pproc_inertial_alfven_wave.py b/examples/DriftKineticMaxwell/inertial_alfven_wave/pproc_inertial_alfven_wave.py new file mode 100644 index 000000000..3e8e0d579 --- /dev/null +++ b/examples/DriftKineticMaxwell/inertial_alfven_wave/pproc_inertial_alfven_wave.py @@ -0,0 +1,95 @@ +"""Post-process the inertial Alfvén wave example: wave frequency from the magnetic energy and energy conservation. + +Run as ``python pproc_inertial_alfven_wave.py [--lam L] [--exb-current]`` (same flags as the params file) or with +the path of the output folder. Writes ``energies.png`` into ``/results`` and prints the fitted frequency. +""" + +import argparse +import sys +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +from scipy.optimize import curve_fit + +from struphy import Output + +SCRIPT_DIR = Path(__file__).resolve().parent +sys.path.insert(0, str(SCRIPT_DIR)) +from params_inertial_alfven_wave import EPSILON, KAPPA, LX, LZ, args # noqa: E402 + +DEFAULT_OUTPUT = Path(args.out_folder) / args.sim_folder + + +def dispersion_roots(lam, exb_current): + """Positive real roots omega of the determinant of the cold-slab dispersion matrix + omega^2 (lam^2 E + chi E_perp) - c1 (k^2 E - k k.E) - P E_z e_z + i omega H E x e_z + (polynomial of degree 6 in omega, interpolated exactly from 7 values).""" + c1, c2, c3 = 1.0, KAPPA**2 * EPSILON, 1.0 / EPSILON + chi, P, H = (EPSILON * KAPPA) ** 2, c2 * c3, c2 if exb_current else 0.0 + k = np.array([2 * np.pi / LX, 0.0, 2 * np.pi / LZ]) + cross_z = np.array([[0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 0.0, 0.0]]) + + def det(w): + mat = ( + w**2 * (lam**2 * np.eye(3) + chi * np.diag([1.0, 1.0, 0.0])) + - c1 * (k @ k * np.eye(3) - np.outer(k, k)) + - P * np.diag([0.0, 0.0, 1.0]) + + 1j * w * H * cross_z + ) + return np.linalg.det(mat).real + + nodes = np.linspace(-1.0, 1.0, 7) + roots = np.roots(np.polyfit(nodes, [det(w) for w in nodes], 6)) + return np.sort(roots[(abs(roots.imag) < 1e-8) & (roots.real > 1e-8)].real) + + +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) + + omega_exact = dispersion_roots(args.lam, args.exb_current)[0] + # en_B ~ cos^2(omega t) = (1 + cos(2 omega t)) / 2 (plus the upper branch for lam > 0, filtered by the fit) + omega = None + if len(t) > 10: + popt, _ = curve_fit( + lambda t, a, b, c, w: a * np.cos(2 * w * t) + b * np.sin(2 * w * t) + c, + t, + en_B / en_B[0], + p0=(0.5, 0.0, 0.5, omega_exact), + ) + omega = abs(popt[3]) + print( + f"fitted omega = {omega:.5f}, analytical omega = {omega_exact:.5f} (rel. error {omega / omega_exact - 1:.2%})" + ) + else: + print("Too few saved time steps for a frequency 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].plot(t, en_B, label="en_B") + axs[0].plot(t, scalars["en_E"].values, label="en_E") + if omega is not None: + axs[0].plot(t, en_B[0] * np.cos(omega_exact * t) ** 2, "r:", label=f"cos^2, omega = {omega_exact:.4f}") + 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"Inertial Alfven wave, {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 3a390184c..a2a51717f 100644 --- a/src/struphy/api/models/__init__.py +++ b/src/struphy/api/models/__init__.py @@ -7,6 +7,7 @@ from struphy.models.drift_kinetic_electrostatic_adiabatic import ( DriftKineticElectrostaticAdiabatic, ) +from struphy.models.drift_kinetic_maxwell import DriftKineticMaxwell from struphy.models.guiding_center import GuidingCenter from struphy.models.gyrokinetic_poisson import GyrokineticPoisson from struphy.models.hasegawa_wakatani import HasegawaWakatani @@ -56,6 +57,7 @@ "CurlCurl", "DeterministicParticleDiffusion", "DriftKineticElectrostaticAdiabatic", + "DriftKineticMaxwell", "GuidingCenter", "GyrokineticPoisson", "HasegawaWakatani", diff --git a/src/struphy/models/__init__.py b/src/struphy/models/__init__.py index ba12868c5..24a3c319f 100644 --- a/src/struphy/models/__init__.py +++ b/src/struphy/models/__init__.py @@ -14,6 +14,7 @@ "CurlCurl": "struphy.models.curl_curl", "DeterministicParticleDiffusion": "struphy.models.deterministic_particle_diffusion", "DriftKineticElectrostaticAdiabatic": "struphy.models.drift_kinetic_electrostatic_adiabatic", + "DriftKineticMaxwell": "struphy.models.drift_kinetic_maxwell", "GuidingCenter": "struphy.models.guiding_center", "GyrokineticPoisson": "struphy.models.gyrokinetic_poisson", "HasegawaWakatani": "struphy.models.hasegawa_wakatani", @@ -60,6 +61,7 @@ from struphy.models.curl_curl import CurlCurl from struphy.models.deterministic_particle_diffusion import DeterministicParticleDiffusion from struphy.models.drift_kinetic_electrostatic_adiabatic import DriftKineticElectrostaticAdiabatic + from struphy.models.drift_kinetic_maxwell import DriftKineticMaxwell from struphy.models.guiding_center import GuidingCenter from struphy.models.gyrokinetic_poisson import GyrokineticPoisson from struphy.models.hasegawa_wakatani import HasegawaWakatani @@ -105,6 +107,7 @@ "CurlCurl", "DeterministicParticleDiffusion", "DriftKineticElectrostaticAdiabatic", + "DriftKineticMaxwell", "GuidingCenter", "GyrokineticPoisson", "HasegawaWakatani", diff --git a/src/struphy/models/drift_kinetic_maxwell.py b/src/struphy/models/drift_kinetic_maxwell.py new file mode 100644 index 000000000..0da63428d --- /dev/null +++ b/src/struphy/models/drift_kinetic_maxwell.py @@ -0,0 +1,345 @@ +import copy +import logging + +from struphy import BaseUnits +from struphy.io.options import LiteralOptions +from struphy.models.base import StruphyModel +from struphy.models.scalars import FunctionScalarFEEC, Scalars +from struphy.models.species import ( + FieldSpecies, + ParticleSpecies, +) +from struphy.models.variables import FEECVariable, PICVariable +from struphy.propagators.drift_kinetic_maxwell_coupling import DriftKineticMaxwellCoupling +from struphy.propagators.push_guiding_center_bx_estar import PushGuidingCenterBxEstar +from struphy.propagators.push_guiding_center_parallel import PushGuidingCenterParallel + +logger = logging.getLogger("struphy") + + +class DriftKineticMaxwell(StruphyModel): + """Electromagnetic drift-kinetic equation for one guiding-center species coupled to Maxwell's equations, in the + thermal normalization, with an implicit, energy conserving field-particle coupling that also holds in the + quasi-neutral limit. + + The electric field enters the guiding-center equations through the gauge-invariant + :math:`\\mathbf E^* = \\mathbf E - \\varepsilon \\mu \\nabla |B|`. Ampère's law contains the guiding-center + current and the (linearized) polarization current; the displacement current has the coefficient + :math:`\\lambda^2 \\varepsilon^2 = \\hat v_\\textrm{th}^2/c^2`, option ``lam`` of the coupling propagator + :class:`~struphy.propagators.drift_kinetic_maxwell_coupling.DriftKineticMaxwellCoupling`; ``lam=0.0`` is the + quasi-neutral limit. + + Parameters + ---------- + base_units: BaseUnits + Base units for normalization (default: ``BaseUnits(kBT=1.0)``; the thermal energy defines the velocity unit). + 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) + 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 + Plasma frequency × length unit / speed of light. If None, computed from units and charge/mass numbers. + """ + + @classmethod + def model_type(cls) -> LiteralOptions.ModelTypes: + return "Kinetic" + + ## species + + class EMFields(FieldSpecies): + def __init__(self): + self.e_field = FEECVariable(space="Hcurl") + self.b_field = FEECVariable(space="Hdiv") + self.init_variables() + + class KineticIons(ParticleSpecies): + def __init__( + self, + charge_number: int = 1, + mass_number: float = 1.0, + epsilon: float = None, + lam: float = None, + kappa: float = None, + ): + self.var = PICVariable(space="Particles5D") + self.init_variables( + charge_number=charge_number, + mass_number=mass_number, + epsilon=epsilon, + lam=lam, + kappa=kappa, + ) + + ## propagators + + class Propagators: + def __init__( + self, + e_field: FEECVariable = None, + b_field: FEECVariable = None, + ): + # the parallel electric force is applied by the coupling step only: push_para gets no phi and no e_field + self.coupling = DriftKineticMaxwellCoupling() + self.push_para = PushGuidingCenterParallel(b_tilde=b_field) + self.push_bxe = PushGuidingCenterBxEstar(b_tilde=b_field, e_field=e_field) + + ## abstract methods + + def __init__( + self, + base_units: BaseUnits = BaseUnits(kBT=1.0), + charge_number: int = 1, + mass_number: float = 1.0, + epsilon: float = None, + lam: float = None, + kappa: float = None, + ): + # 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, 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(e_field=self.em_fields.e_field, b_field=self.em_fields.b_field) + + # 4. assign variables to propagators + 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_para.variables.ions = self.kinetic_ions.var + self.propagators.push_bxe.variables.ions = self.kinetic_ions.var + + # 5. define scalars to be tracked during simulation: the parts of the energy conserved by the coupling step + # (with the coefficients 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]) + self.scalars = Scalars( + en_E=electric_energy, + en_B=magnetic_energy, + en_f=particle_energy, + en_tot=electric_energy + magnetic_energy + particle_energy, + ) + + @property + def bulk_species(self): + return self.kinetic_ions + + @property + def velocity_scale(self): + return "thermal" + + def post_allocate(self): + """Check that the parallel electric force is applied only once (by the coupling step), then no initial field + solve: the electric field is initialized from the parameter file. + + :meta private: + """ + push_para = self.propagators.push_para + assert push_para.e_field is None and push_para._allocate_phi, ( + "DriftKineticMaxwell: the parallel electric force b* . E is applied by the coupling step; " + "push_para must not get an electric field or a potential (it would be applied twice)." + ) + + ## default parameters + + def generate_default_parameter_file(self, path=None, prompt=True): + """The thermal velocity scale needs the thermal energy unit: ``BaseUnits(kBT=1.0)``. + + :meta private: + """ + params_path = super().generate_default_parameter_file(path=path, prompt=prompt) + with open(params_path, "r") as f: + lines = f.readlines() + with open(params_path, "w") as f: + for line in lines: + f.write("base_units = BaseUnits(kBT=1.0)\n" if "BaseUnits(" in line else line) + return params_path + + @classmethod + def doc_pde(cls): + r"""**PDEs solved by model:** + + Drift-kinetic equation for :math:`f(\mathbf X, v_\parallel, \mu, t)`: + + .. math:: + + \frac{\partial f}{\partial t} + \frac{v_\parallel \mathbf B^* + \mathbf E^* \times \mathbf b_0}{B^*_\parallel} + \cdot \nabla f + \frac{1}{\varepsilon} \frac{\mathbf B^* \cdot \mathbf E^*}{B^*_\parallel} + \frac{\partial f}{\partial v_\parallel} = 0\,, + + with + + .. math:: + + \mathbf E^* = \mathbf E - \varepsilon \mu \nabla |\mathbf B|\,,\qquad + \mathbf B^* = \mathbf B_0 + \varepsilon v_\parallel \nabla \times \mathbf b_0\,,\qquad + B^*_\parallel = \mathbf b_0 \cdot (\mathbf B_0 + \tilde{\mathbf B}) + \varepsilon v_\parallel + (\nabla \times \mathbf b_0) \cdot \mathbf b_0\,, + + where :math:`|\mathbf B| = \mathbf b_0 \cdot (\mathbf B_0 + \tilde{\mathbf B})`. + + Ampère's law, with the displacement current and the (linearized) polarization current of the guiding centers: + + .. math:: + + \varepsilon^2 \lambda^2 \frac{\partial \mathbf E}{\partial t} + + \varepsilon^2 \kappa^2 \frac{n_0}{|\mathbf B_0|^2} \frac{\partial \mathbf E_\perp}{\partial t} + = \nabla \times \tilde{\mathbf B} - \varepsilon \kappa^2 \iint f B^*_\parallel + \left( v_\parallel \frac{\mathbf B^*}{B^*_\parallel} + \frac{\mathbf E \times \mathbf b_0}{B^*_\parallel} + \right) \textrm d v_\parallel \textrm d \mu\,, + + Faraday's law: + + .. math:: + + \frac{\partial \tilde{\mathbf B}}{\partial t} + \nabla \times \mathbf E = 0\,. + + :math:`\mathbf B_0` is a static equilibrium field, :math:`n_0` the equilibrium density, and + :math:`\mathbf E_\perp = \mathbf E - (\mathbf b_0 \cdot \mathbf E) \mathbf b_0`. The current in Ampère's law is + the guiding-center current without the magnetic (grad-:math:`B`) drift, whose work on :math:`\mathbf E` is + therefore missing where :math:`\nabla |B| \neq 0`. The coefficients :math:`\varepsilon\lambda` (option ``lam``, + 0 is the quasi-neutral limit), :math:`(\varepsilon\kappa)^2` (option ``chi``) and :math:`\varepsilon\kappa^2` + (option ``c2``) can be changed in + :class:`~struphy.propagators.drift_kinetic_maxwell_coupling.DriftKineticMaxwellCoupling`. Without the + polarization current the quasi-neutral limit would be singular for flute modes (:math:`k_\parallel = 0`). + """ + + @classmethod + def doc_normalization(cls): + r"""The model uses the thermal velocity of the species as reference velocity: + + .. math:: + + \hat v = \hat v_\textrm{th} = \sqrt{\frac{k_\textrm{B} \hat T}{m}}\,,\qquad \hat t = \frac{\hat x}{\hat v}\,, + \qquad \hat E = \hat v \hat B\,,\qquad \hat \mu = \frac{m \hat v^2}{\hat B}\,. + + The species parameters (see :ref:`equation_params`) are + + .. math:: + + \varepsilon = \frac{1}{\hat\Omega_c \hat t} = \frac{\hat v}{\hat\Omega_c \hat x}\,,\qquad + \lambda = \frac{\hat\Omega_c \hat x}{c}\,,\qquad \kappa = \frac{\hat\Omega_p \hat x}{c}\,. + + The drift-kinetic equation and Faraday's law then have only the coefficient :math:`\varepsilon` + (:math:`\frac{m \hat v}{q \hat B \hat x} = \varepsilon` in :math:`\mathbf B^*` and :math:`\mathbf E^*`, + :math:`\frac{q \hat B \hat t}{m} = \frac{1}{\varepsilon}` in the parallel force). Ampère's law + :math:`\frac{1}{c^2} \partial_t \mathbf E + \mu_0 \frac{n m}{|B|^2} \partial_t \mathbf E_\perp = \nabla \times + \mathbf B - \mu_0 q \mathbf J`, multiplied by :math:`\hat x / \hat B`, has the coefficients + + .. math:: + + \frac{\hat v^2}{c^2} = \varepsilon^2 \lambda^2\,,\qquad + \frac{\mu_0 \hat n m \hat v^2}{\hat B^2} = \frac{\hat v^2}{\hat v_A^2} = \varepsilon^2 \kappa^2\,,\qquad + \frac{\mu_0 q \hat n \hat v \hat x}{\hat B} = \varepsilon \kappa^2\,. + + Their ratio :math:`\kappa^2/\lambda^2 = \hat\Omega_p^2/\hat\Omega_c^2 = c^2/\hat v_A^2` compares the polarization + with the vacuum dielectric, and :math:`1/\kappa` is the inertial length :math:`c/\hat\Omega_p` in units of + :math:`\hat x`. For example, deuterium with :math:`\hat x = 1` m, :math:`\hat B = 5` T, + :math:`\hat n = 10^{20}` m\ :sup:`-3`, :math:`k_\textrm{B} \hat T = 10` keV gives :math:`\varepsilon = 2.9 \cdot + 10^{-3}`, :math:`\lambda = 0.80`, :math:`\kappa = 31`, i.e. :math:`\varepsilon^2 \lambda^2 = 5.3 \cdot 10^{-6}`, + :math:`\varepsilon^2 \kappa^2 = 8.0 \cdot 10^{-3}` and :math:`\varepsilon \kappa^2 = 2.8`.""" + + @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 of + :class:`~struphy.propagators.drift_kinetic_maxwell_coupling.DriftKineticMaxwellCoupling`, here the defaults): + + - Electric field energy (displacement and polarization) + :math:`\frac{\varepsilon^2 \lambda^2}{2} \| \mathbf E \|^2 + \frac{\varepsilon^2 \kappa^2}{2} + \| \mathbf E_\perp \|^2_{n_0/|B_0|^2}`: ``en_E`` + - Magnetic field energy :math:`\frac{1}{2} \| \tilde{\mathbf B} \|^2`: ``en_B`` + - Particle energy :math:`\varepsilon^2 \kappa^2 \sum_p w_p (v_{\parallel,p}^2/2 + \mu_p |B_0|)`: ``en_f`` + - Total energy: ``en_tot`` (conserved by the coupling step up to the solver tolerance)""" + + @classmethod + def doc_discretization(cls): + """Time integration is performed by the following propagators (in sequence): + + 1. :class:`~struphy.propagators.drift_kinetic_maxwell_coupling.DriftKineticMaxwellCoupling` + 2. :class:`~struphy.propagators.push_guiding_center_parallel.PushGuidingCenterParallel` + 3. :class:`~struphy.propagators.push_guiding_center_bx_estar.PushGuidingCenterBxEstar` + """ + doc = rf"""Each term of the drift-kinetic equation is applied by exactly one propagator; the parallel electric +force only by the coupling step (the model passes no potential and no electric field to ``push_para``, and asserts +this): + +- :math:`\dot v_\parallel \ni \frac{{1}}{{\varepsilon}} \mathbf b^* \cdot \mathbf E`: coupling step + (``push_vpar_with_efield``); +- :math:`\dot v_\parallel \ni -\mu\, \mathbf b^* \cdot \nabla |B|` (mirror force): ``push_para``; +- :math:`\dot{{\mathbf X}} \ni v_\parallel \mathbf b^*` (streaming, curvature drift): ``push_para``; its current + enters Ampère's law in the coupling step; +- :math:`\dot{{\mathbf X}} \ni \mathbf E \times \mathbf b_0 / B^*_\parallel`: ``push_bxe``; its current enters + Ampère's law in the coupling step (option ``exb_current``); +- :math:`\dot{{\mathbf X}} \ni -\varepsilon \mu \nabla |B| \times \mathbf b_0 / B^*_\parallel` (grad-:math:`B` + drift): ``push_bxe``; its current is not included. + +**1. DriftKineticMaxwellCoupling:** + +{DriftKineticMaxwellCoupling.__doc__} + +**2. PushGuidingCenterParallel:** + +{PushGuidingCenterParallel.__doc__} + +**3. PushGuidingCenterBxEstar:** + +{PushGuidingCenterBxEstar.__doc__} +""" + return doc + + @classmethod + def doc_long_description(cls): + r"""DriftKineticMaxwell couples one drift-kinetic (guiding-center) species to Maxwell's equations in the + gauge-invariant :math:`v_\parallel` formulation (the electric field, not a potential, enters + :math:`\mathbf E^*`), normalized with the thermal velocity of the species. Maxwell's equations and the parallel + acceleration are solved together in one implicit, energy conserving step, which allows time steps above the + plasma period and has no cancellation problem. Since the drift-kinetic ordering removes the polarization drift, + its (linearized) polarization current is added to Ampère's law; with it, the quasi-neutral limit without + displacement current is well-posed and gives the inertial Alfvén wave. The guiding centers then stream along + the field and drift with :math:`\mathbf E^* \times \mathbf b_0`. The total energy is conserved exactly by the + splitting if :math:`\nabla |B_0| = 0` and either :math:`\mu = 0` or :math:`\mathbf b_0 \cdot \tilde{\mathbf B} + = 0`.""" + + @classmethod + def doc_examples(cls): + r"""Create a drift-kinetic Maxwell model in the quasi-neutral limit: + + .. code-block:: python + + from struphy.models import DriftKineticMaxwell + + model = DriftKineticMaxwell() + prop = model.propagators.coupling + prop.options = prop.Options(lam=0.0) + """ + + @classmethod + def doc_use_cases(cls): + r"""This model is appropriate for: + + - electromagnetic drift-kinetic simulations with time steps above the plasma period + - shear (inertial) Alfvén waves with kinetic parallel dynamics + - the quasi-neutral limit (``lam=0.0``) with full-f weights""" + + @classmethod + def doc_cannot_be_used_for(cls): + r"""This model is not suitable for: + + - polar splines (not supported yet by the coupling propagator) + - the quasi-neutral limit with control-variate (delta-f) weights or with vacuum regions + - finite Larmor radius effects beyond the linearized long-wavelength polarization + - exact energy conservation in inhomogeneous :math:`|B_0|` (the grad-:math:`B` drift current is missing)""" diff --git a/src/struphy/models/tests/verification/test_verif_DriftKineticMaxwell.py b/src/struphy/models/tests/verification/test_verif_DriftKineticMaxwell.py new file mode 100644 index 000000000..a33dcf944 --- /dev/null +++ b/src/struphy/models/tests/verification/test_verif_DriftKineticMaxwell.py @@ -0,0 +1,144 @@ +"""Verification of DriftKineticMaxwell: the shear (inertial) Alfvén wave of a cold guiding-center plasma in a slab +with uniform B0 = e_z, against the roots of the dispersion relation + +.. math:: + + \\omega^2 (\\lambda^2 \\mathbf E + \\chi \\mathbf E_\\perp) = c_1 (k^2 \\mathbf E - \\mathbf k\\, \\mathbf k \\cdot \\mathbf E) + + P E_z \\mathbf e_z - \\mathrm i \\omega H\\, \\mathbf E \\times \\mathbf e_z\\,, + +with :math:`P = c_2 c_3 n_0` (parallel response) and :math:`H = c_2 n_0` (E x B current, zero if switched off), in the +thermal normalization: :math:`\\chi = (\\varepsilon\\kappa_\\textrm{s})^2`, :math:`c_2 = \\varepsilon\\kappa_\\textrm{s}^2`, +:math:`c_3 = 1/\\varepsilon`, and :math:`\\lambda` the option ``lam`` (coefficient of the displacement current). +For :math:`\\lambda = 0` and :math:`H = 0` the shear branch is the inertial Alfvén wave +:math:`\\omega^2 = (c_1/\\chi)\\, k_z^2 / (1 + k_x^2 c_1/P)`. The energy is conserved exactly by the splitting +(cold plasma, uniform B0). + +Run it on several processes with:: + + mpirun -n 4 pytest -p no:xdist --with-mpi src/struphy/models/tests/verification/test_verif_DriftKineticMaxwell.py +""" + +import logging + +import numpy as np +import pytest + +from struphy import set_logging_level + +set_logging_level(logging.WARNING) + +# species parameters (in the thermal normalization of the model): chi = (EPSILON KAPPA)^2 = 9, c2 = 18, c3 = 2 +EPSILON, LAM, KAPPA = 0.5, 2.0, 6.0 +LX, LZ = 1.0, 4.0 +VTH = 1e-6 # cold + + +def _dispersion_roots(lam, exb_current): + """Positive real roots omega of the determinant of the cold-slab dispersion matrix (polynomial of degree 6, + interpolated exactly from 7 values).""" + c1, c2, c3 = 1.0, KAPPA**2 * EPSILON, 1.0 / EPSILON + chi, P, H = (EPSILON * KAPPA) ** 2, c2 * c3, c2 if exb_current else 0.0 + k = np.array([2 * np.pi / LX, 0.0, 2 * np.pi / LZ]) + proj_perp = np.diag([1.0, 1.0, 0.0]) + cross_z = np.array([[0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 0.0, 0.0]]) # E -> E x e_z + + def det(w): + mat = ( + w**2 * (lam**2 * np.eye(3) + chi * proj_perp) + - c1 * (k @ k * np.eye(3) - np.outer(k, k)) + - P * np.diag([0.0, 0.0, 1.0]) + + 1j * w * H * cross_z + ) + return np.linalg.det(mat) + + nodes = np.linspace(-1.0, 1.0, 7) + coeffs = np.polyfit(nodes, [det(w).real for w in nodes], 6) + roots = np.roots(coeffs) + return np.sort(roots[(abs(roots.imag) < 1e-8) & (roots.real > 1e-8)].real) + + +def _alfven_signal(lam, exb_current, dt, n_steps): + """Time series s(t) = b(t)^T M2 b(0) / |b(0)|^2 of the shear Alfvén wave, excited by a B_y perturbation, and the + relative change of the energy.""" + from struphy import ( + DerhamOptions, + LoadingParameters, + Simulation, + domains, + equils, + grids, + maxwellians, + models, + perturbations, + ) + from struphy.linear_algebra.solver import SolverParameters + + model = models.DriftKineticMaxwell(epsilon=EPSILON, lam=LAM, kappa=KAPPA) + model.kinetic_ions.set_markers( + loading_params=LoadingParameters(ppc=50, seed=1234, moments=(0.0, 0.0, VTH, VTH), B0=1.0) + ) + cold = dict(n=(1.0, None), vth_para=(VTH, None), vth_perp=(VTH, None), B0=1.0) + model.kinetic_ions.var.add_background(maxwellians.GyroMaxwellian2D(**cold)) + model.kinetic_ions.var.add_initial_condition(maxwellians.GyroMaxwellian2D(**cold)) + model.em_fields.b_field.add_perturbation( + perturbations.ModesSin(ls=(1,), ns=(1,), given_in_basis="2", comp=1, amps=(1e-2,)) + ) + prop = model.propagators.coupling + solver = {} if exb_current else {"solver": "pcg"} + prop.options = prop.Options( + lam=lam, exb_current=exb_current, solver_params=SolverParameters(tol=1e-12, maxiter=3000), **solver + ) + Simulation( + model=model, + domain=domains.Cuboid(r1=LX, r2=1.0, r3=LZ), + equil=equils.HomogenSlab(B0z=1.0), + grid=grids.TensorProductGrid(num_elements=(8, 1, 8)), + derham_opts=DerhamOptions(degree=(3, 1, 3)), + ).allocate() + + props = model.propagators + prop.prepare_time_step(dt) # the initial projection for lam = 0 + energy_0 = prop.energy() + b = model.em_fields.b_field.spline.vector + b_init = b.copy() + m2_b_init = prop._M2.dot(b_init) + signal = [b.inner(m2_b_init)] + for _ in range(n_steps): + props.coupling(dt) + props.push_para(dt) + props.push_bxe(dt) + signal.append(b.inner(m2_b_init)) + return dt * np.arange(n_steps + 1), np.array(signal) / signal[0], prop.energy() / energy_0 - 1.0 + + +def _fit(t, s, omegas): + """Least-squares fit of s(t) by a constant plus one harmonic per initial guess in omegas.""" + from scipy.optimize import curve_fit + + def model(t, c, *p): + out = c + 0.0 * t + for i in range(len(p) // 3): + out = out + p[3 * i] * np.cos(p[3 * i + 2] * t) + p[3 * i + 1] * np.sin(p[3 * i + 2] * t) + return out + + p0 = [0.0] + for w in omegas: + p0 += [0.5, 0.0, w] + popt, _ = curve_fit(model, t, s, p0=p0) + return np.abs(popt[3::3]) + + +@pytest.mark.parametrize( + "lam, exb_current, dt, n_steps", + [(0.0, False, 0.5, 160), (0.0, True, 0.5, 160), (1.0, False, 0.1, 400)], +) +def test_shear_alfven_wave(lam, exb_current, dt, n_steps): + roots = _dispersion_roots(lam, exb_current) + t, s, energy_change = _alfven_signal(lam, exb_current, dt, n_steps) + if lam > 0.0: + # the B_y perturbation also excites the fast (light-wave like) upper branch of the shear mode + omegas = _fit(t, s, (roots[0], roots[-1])) + else: + omegas = _fit(t, s, (roots[0],)) + assert abs(omegas[0] / roots[0] - 1.0) < 0.02, (omegas, roots) + assert abs(energy_change) < 1e-7 diff --git a/src/struphy/pic/accumulation/kernels/dk_maxwell/__init__.py b/src/struphy/pic/accumulation/kernels/dk_maxwell/__init__.py new file mode 100644 index 000000000..9039bc198 --- /dev/null +++ b/src/struphy/pic/accumulation/kernels/dk_maxwell/__init__.py @@ -0,0 +1,10 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported.""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +# the arguments the kernel writes to, by position: mat11 to vec3 +OUTPUTS = tuple(range(3, 15)) + +dk_maxwell = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, host_options={"outputs": OUTPUTS}) diff --git a/src/struphy/pic/accumulation/kernels/dk_maxwell/dk_maxwell_kernels.py b/src/struphy/pic/accumulation/kernels/dk_maxwell/dk_maxwell_kernels.py new file mode 100644 index 000000000..9b1eee3a9 --- /dev/null +++ b/src/struphy/pic/accumulation/kernels/dk_maxwell/dk_maxwell_kernels.py @@ -0,0 +1,156 @@ +"Accumulation kernel for gyro-center (5D) particles." + +from numpy import empty +from pyccel.decorators import stack_array + +import struphy.geometry.evaluation_kernels as evaluation_kernels +import struphy.kernel_arguments.pusher_args_kernels as pusher_args_kernels # do not remove; needed to identify dependencies +import struphy.linear_algebra.linalg_kernels as linalg_kernels +import struphy.pic.accumulation.particle_to_mat_kernels as particle_to_mat_kernels +from struphy.bsplines.evaluation_kernels_3d import ( + eval_0form_spline_mpi, + eval_1form_spline_mpi, + eval_2form_spline_mpi, + get_spans, +) +from struphy.kernel_arguments.pusher_args_kernels import DerhamArguments, DomainArguments, MarkerArguments + + +@stack_array("dfm", "b2", "b_star", "unit_b1", "a", "filling_m") +def dk_maxwell( + args_markers: "MarkerArguments", + args_derham: "DerhamArguments", + args_domain: "DomainArguments", + mat11: "float[:,:,:,:,:,:]", + mat12: "float[:,:,:,:,:,:]", + mat13: "float[:,:,:,:,:,:]", + mat21: "float[:,:,:,:,:,:]", + mat22: "float[:,:,:,:,:,:]", + mat23: "float[:,:,:,:,:,:]", + mat31: "float[:,:,:,:,:,:]", + mat32: "float[:,:,:,:,:,:]", + mat33: "float[:,:,:,:,:,:]", + vec1: "float[:,:,:]", + vec2: "float[:,:,:]", + vec3: "float[:,:,:]", + epsilon: float, + skew: float, + b2_1: "float[:,:,:]", + b2_2: "float[:,:,:]", + b2_3: "float[:,:,:]", + curl_unit_b2_1: "float[:,:,:]", + curl_unit_b2_2: "float[:,:,:]", + curl_unit_b2_3: "float[:,:,:]", + unit_b1_1: "float[:,:,:]", + unit_b1_2: "float[:,:,:]", + unit_b1_3: "float[:,:,:]", + B_dot_b_coeffs: "float[:,:,:]", + curl_unit_b_dot_b0: "float[:,:,:]", +): + r"""Accumulation kernel for the propagator + :class:`~struphy.propagators.drift_kinetic_maxwell_coupling.DriftKineticMaxwellCoupling`. Accumulates into + :math:`V^1_h` with the (non-symmetric) filling functions + + .. math:: + + A_p^{\mu\nu} = w_p \left( a_{p,\mu} a_{p,\nu} + \sigma\, S_{p,\mu\nu} \right)\,, + \qquad + B_p^\mu = w_p\, v_{\parallel,p}\, a_{p,\mu}\,, + \qquad + \mathbf a_p = \frac{\hat{\mathbf B}^{*2}}{\hat B^{*3}_\parallel}(\boldsymbol \eta_p)\,, + \qquad + S_p \mathbf x = \frac{\mathbf x \times \hat{\mathbf b}^1_0}{\hat B^{*3}_\parallel}(\boldsymbol \eta_p)\,, + + with :math:`\hat{\mathbf B}^{*2} = \hat{\mathbf B}^2 + \varepsilon v_\parallel \hat \nabla \times + \hat{\mathbf b}^1_0` and :math:`\hat B^{*3}_\parallel = \sqrt g\, (\hat B_\parallel + \varepsilon v_\parallel + \widehat{(\nabla \times \mathbf b_0) \cdot \mathbf b_0})`, as in + :class:`~struphy.propagators.push_guiding_center_parallel.PushGuidingCenterParallel` and + :class:`~struphy.propagators.push_guiding_center_bx_estar.PushGuidingCenterBxEstar`. The rank-one part is the + parallel response, the antisymmetric part :math:`S_p` the :math:`E \times B` current, scaled by ``skew`` + (:math:`\sigma`). + """ + + markers = args_markers.markers + + dfm = empty((3, 3), dtype=float) + b2 = empty(3, dtype=float) + b_star = empty(3, dtype=float) + unit_b1 = empty(3, dtype=float) + a = empty(3, dtype=float) + filling_m = empty((3, 3), dtype=float) + + for ip in range(args_markers.n_markers): + # only do something if particle is a "true" particle (i.e. not a hole) + if markers[0, ip] == -1.0: + continue + + eta1 = markers[0, ip] + eta2 = markers[1, ip] + eta3 = markers[2, ip] + v = markers[3, ip] + weight = markers[5, ip] + + evaluation_kernels.df(eta1, eta2, eta3, args_domain, dfm) + det_df = linalg_kernels.det(dfm) + + span1, span2, span3 = get_spans(eta1, eta2, eta3, args_derham) + + # B* (2-form) and B*_parallel (3-form) + eval_2form_spline_mpi(span1, span2, span3, args_derham, b2_1, b2_2, b2_3, b2) + eval_2form_spline_mpi(span1, span2, span3, args_derham, curl_unit_b2_1, curl_unit_b2_2, curl_unit_b2_3, b_star) + b_star *= epsilon * v + b_star += b2 + + b_star_parallel = eval_0form_spline_mpi(span1, span2, span3, args_derham, curl_unit_b_dot_b0) + b_star_parallel *= epsilon * v + b_star_parallel += eval_0form_spline_mpi(span1, span2, span3, args_derham, B_dot_b_coeffs) + b_star_parallel *= det_df + + a[:] = b_star / b_star_parallel + + # unit vector b0 (1-form) + eval_1form_spline_mpi(span1, span2, span3, args_derham, unit_b1_1, unit_b1_2, unit_b1_3, unit_b1) + + # filling: w (a a^T + skew S), with S x = x cross b0 / B*_parallel + for i in range(3): + for j in range(3): + filling_m[i, j] = weight * a[i] * a[j] + + s = weight * skew / b_star_parallel + filling_m[0, 1] += s * unit_b1[2] + filling_m[0, 2] -= s * unit_b1[1] + filling_m[1, 0] -= s * unit_b1[2] + filling_m[1, 2] += s * unit_b1[0] + filling_m[2, 0] += s * unit_b1[1] + filling_m[2, 1] -= s * unit_b1[0] + + particle_to_mat_kernels.m_v_fill_b_v1_full( + args_derham, + eta1, + eta2, + eta3, + mat11, + mat12, + mat13, + mat21, + mat22, + mat23, + mat31, + mat32, + mat33, + filling_m[0, 0], + filling_m[0, 1], + filling_m[0, 2], + filling_m[1, 0], + filling_m[1, 1], + filling_m[1, 2], + filling_m[2, 0], + filling_m[2, 1], + filling_m[2, 2], + vec1, + vec2, + vec3, + weight * v * a[0], + weight * v * a[1], + weight * v * a[2], + ) diff --git a/src/struphy/pic/accumulation/kernels/dk_maxwell_mf/__init__.py b/src/struphy/pic/accumulation/kernels/dk_maxwell_mf/__init__.py new file mode 100644 index 000000000..f1daa91ad --- /dev/null +++ b/src/struphy/pic/accumulation/kernels/dk_maxwell_mf/__init__.py @@ -0,0 +1,10 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported.""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +# the arguments the kernel writes to, by position: args_markers, vec1, vec2, vec3 +OUTPUTS = (0, 3, 4, 5) + +dk_maxwell_mf = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, host_options={"outputs": OUTPUTS}) diff --git a/src/struphy/pic/accumulation/kernels/dk_maxwell_mf/dk_maxwell_mf_kernels.py b/src/struphy/pic/accumulation/kernels/dk_maxwell_mf/dk_maxwell_mf_kernels.py new file mode 100644 index 000000000..f869e9e42 --- /dev/null +++ b/src/struphy/pic/accumulation/kernels/dk_maxwell_mf/dk_maxwell_mf_kernels.py @@ -0,0 +1,117 @@ +"Accumulation kernel for gyro-center (5D) particles, matrix-free variant of dk_maxwell." + +from numpy import empty +from pyccel.decorators import stack_array + +import struphy.geometry.evaluation_kernels as evaluation_kernels +import struphy.kernel_arguments.pusher_args_kernels as pusher_args_kernels # do not remove; needed to identify dependencies +import struphy.linear_algebra.linalg_kernels as linalg_kernels +import struphy.pic.accumulation.particle_to_mat_kernels as particle_to_mat_kernels +from struphy.bsplines.evaluation_kernels_3d import ( + eval_0form_spline_mpi, + eval_1form_spline_mpi, + eval_2form_spline_mpi, + get_spans, +) +from struphy.kernel_arguments.pusher_args_kernels import DerhamArguments, DomainArguments, MarkerArguments + + +@stack_array("dfm", "b2", "b_star", "unit_b1", "a") +def dk_maxwell_mf( + args_markers: "MarkerArguments", + args_derham: "DerhamArguments", + args_domain: "DomainArguments", + vec1: "float[:,:,:]", + vec2: "float[:,:,:]", + vec3: "float[:,:,:]", + epsilon: float, + b2_1: "float[:,:,:]", + b2_2: "float[:,:,:]", + b2_3: "float[:,:,:]", + curl_unit_b2_1: "float[:,:,:]", + curl_unit_b2_2: "float[:,:,:]", + curl_unit_b2_3: "float[:,:,:]", + unit_b1_1: "float[:,:,:]", + unit_b1_2: "float[:,:,:]", + unit_b1_3: "float[:,:,:]", + B_dot_b_coeffs: "float[:,:,:]", + curl_unit_b_dot_b0: "float[:,:,:]", + column: "int", +): + r"""Accumulation kernel for the propagator + :class:`~struphy.propagators.drift_kinetic_maxwell_coupling.DriftKineticMaxwellCoupling` (matrix-free). + Accumulates the vector :math:`B_p^\mu = w_p v_{\parallel,p} a_{p,\mu}` into :math:`V^1_h` and stores, per marker, + the data of the two particle matrices of + :func:`~struphy.pic.accumulation.kernels.dk_maxwell.dk_maxwell_kernels.dk_maxwell` + (see :class:`~struphy.pic.accumulation.particles_to_grid.MarkerOperator`): + + - the rank-one filling :math:`w_p\, \mathbf a_p \mathbf a_p^\top` in the columns ``column`` to ``column + 3`` + (:math:`\mathbf a_p` and :math:`w_p`), applied by + :func:`~struphy.pic.accumulation.kernels.marker_operator_rank1_dot.marker_operator_rank1_dot_kernels.marker_operator_rank1_dot`; + - the antisymmetric filling :math:`w_p S_p`, :math:`S_p \mathbf x = \mathbf x \times \hat{\mathbf b}^1_0 / + \hat B^{*3}_\parallel`, entries 12, 13, 23 in the columns ``column + 4`` to ``column + 6``, applied by + :func:`~struphy.pic.accumulation.kernels.marker_operator_asym_dot.marker_operator_asym_dot_kernels.marker_operator_asym_dot`. + """ + + markers = args_markers.markers + + dfm = empty((3, 3), dtype=float) + b2 = empty(3, dtype=float) + b_star = empty(3, dtype=float) + unit_b1 = empty(3, dtype=float) + a = empty(3, dtype=float) + + for ip in range(args_markers.n_markers): + # only do something if particle is a "true" particle (i.e. not a hole) + if markers[0, ip] == -1.0: + continue + + eta1 = markers[0, ip] + eta2 = markers[1, ip] + eta3 = markers[2, ip] + v = markers[3, ip] + weight = markers[5, ip] + + evaluation_kernels.df(eta1, eta2, eta3, args_domain, dfm) + det_df = linalg_kernels.det(dfm) + + span1, span2, span3 = get_spans(eta1, eta2, eta3, args_derham) + + # B* (2-form) and B*_parallel (3-form) + eval_2form_spline_mpi(span1, span2, span3, args_derham, b2_1, b2_2, b2_3, b2) + eval_2form_spline_mpi(span1, span2, span3, args_derham, curl_unit_b2_1, curl_unit_b2_2, curl_unit_b2_3, b_star) + b_star *= epsilon * v + b_star += b2 + + b_star_parallel = eval_0form_spline_mpi(span1, span2, span3, args_derham, curl_unit_b_dot_b0) + b_star_parallel *= epsilon * v + b_star_parallel += eval_0form_spline_mpi(span1, span2, span3, args_derham, B_dot_b_coeffs) + b_star_parallel *= det_df + + a[:] = b_star / b_star_parallel + + # unit vector b0 (1-form) + eval_1form_spline_mpi(span1, span2, span3, args_derham, unit_b1_1, unit_b1_2, unit_b1_3, unit_b1) + + s = weight / b_star_parallel + + markers[column, ip] = a[0] + markers[column + 1, ip] = a[1] + markers[column + 2, ip] = a[2] + markers[column + 3, ip] = weight + markers[column + 4, ip] = s * unit_b1[2] + markers[column + 5, ip] = -s * unit_b1[1] + markers[column + 6, ip] = s * unit_b1[0] + + particle_to_mat_kernels.vec_fill_v1( + args_derham, + span1, + span2, + span3, + vec1, + vec2, + vec3, + weight * v * a[0], + weight * v * a[1], + weight * v * a[2], + ) diff --git a/src/struphy/pic/pushing/kernels/push_vpar_with_efield/__init__.py b/src/struphy/pic/pushing/kernels/push_vpar_with_efield/__init__.py new file mode 100644 index 000000000..e6dd0abaf --- /dev/null +++ b/src/struphy/pic/pushing/kernels/push_vpar_with_efield/__init__.py @@ -0,0 +1,10 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported.""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +# the arguments the kernel writes to, by position: args_markers +OUTPUTS = (2,) + +push_vpar_with_efield = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, host_options={"outputs": OUTPUTS}) diff --git a/src/struphy/pic/pushing/kernels/push_vpar_with_efield/push_vpar_with_efield_kernels.py b/src/struphy/pic/pushing/kernels/push_vpar_with_efield/push_vpar_with_efield_kernels.py new file mode 100644 index 000000000..e5a5e7b11 --- /dev/null +++ b/src/struphy/pic/pushing/kernels/push_vpar_with_efield/push_vpar_with_efield_kernels.py @@ -0,0 +1,90 @@ +"Pusher kernel for gyro-center (5D) particles." + +from numpy import empty +from pyccel.decorators import stack_array + +import struphy.geometry.evaluation_kernels as evaluation_kernels +import struphy.kernel_arguments.pusher_args_kernels as pusher_args_kernels # do not remove; needed to identify dependencies +import struphy.linear_algebra.linalg_kernels as linalg_kernels +from struphy.bsplines.evaluation_kernels_3d import ( + eval_0form_spline_mpi, + eval_1form_spline_mpi, + eval_2form_spline_mpi, + get_spans, +) +from struphy.kernel_arguments.pusher_args_kernels import DerhamArguments, DomainArguments, MarkerArguments + + +@stack_array("dfm", "b2", "b_star", "e_form") +def push_vpar_with_efield( + dt: float, + stage: int, + args_markers: "MarkerArguments", + args_domain: "DomainArguments", + args_derham: "DerhamArguments", + epsilon: float, + e1_1: "float[:,:,:]", + e1_2: "float[:,:,:]", + e1_3: "float[:,:,:]", + b2_1: "float[:,:,:]", + b2_2: "float[:,:,:]", + b2_3: "float[:,:,:]", + curl_unit_b2_1: "float[:,:,:]", + curl_unit_b2_2: "float[:,:,:]", + curl_unit_b2_3: "float[:,:,:]", + B_dot_b_coeffs: "float[:,:,:]", + curl_unit_b_dot_b0: "float[:,:,:]", + const_factor: float, +): + r"""Updates the parallel velocities of guiding centers at fixed positions as + + .. math:: + + \frac{v_{\parallel,p}^{n+1} - v_{\parallel,p}^n}{\Delta t} + = c\, \frac{\hat{\mathbf B}^{*2}}{\hat B^{*3}_\parallel} \cdot \hat{\mathbf E}^1 (\boldsymbol \eta_p)\,, + + where :math:`\hat{\mathbf E}^1` is the 1-form spline field of the given FE coefficients, :math:`c` + (``const_factor``) a constant, and :math:`\hat{\mathbf B}^{*2}`, :math:`\hat B^{*3}_\parallel` are evaluated with + :math:`v_{\parallel,p}^n` as in + :func:`~struphy.pic.accumulation.kernels.dk_maxwell.dk_maxwell_kernels.dk_maxwell`. + """ + + dfm = empty((3, 3), dtype=float) + b2 = empty(3, dtype=float) + b_star = empty(3, dtype=float) + e_form = empty(3, dtype=float) + + markers = args_markers.markers + n_markers = args_markers.n_markers + valid_mks = args_markers.valid_mks + + for ip in range(n_markers): + # only do something if particle is valid (i.e. not a hole or ghost) + if not valid_mks[ip]: + continue + + eta1 = markers[0, ip] + eta2 = markers[1, ip] + eta3 = markers[2, ip] + v = markers[3, ip] + + evaluation_kernels.df(eta1, eta2, eta3, args_domain, dfm) + det_df = linalg_kernels.det(dfm) + + span1, span2, span3 = get_spans(eta1, eta2, eta3, args_derham) + + # B* (2-form) and B*_parallel (3-form) + eval_2form_spline_mpi(span1, span2, span3, args_derham, b2_1, b2_2, b2_3, b2) + eval_2form_spline_mpi(span1, span2, span3, args_derham, curl_unit_b2_1, curl_unit_b2_2, curl_unit_b2_3, b_star) + b_star *= epsilon * v + b_star += b2 + + b_star_parallel = eval_0form_spline_mpi(span1, span2, span3, args_derham, curl_unit_b_dot_b0) + b_star_parallel *= epsilon * v + b_star_parallel += eval_0form_spline_mpi(span1, span2, span3, args_derham, B_dot_b_coeffs) + b_star_parallel *= det_df + + # electric field (1-form) + eval_1form_spline_mpi(span1, span2, span3, args_derham, e1_1, e1_2, e1_3, e_form) + + markers[3, ip] += dt * const_factor * linalg_kernels.scalar_dot(b_star, e_form) / b_star_parallel diff --git a/src/struphy/pic/tests/test_kernel_backends.py b/src/struphy/pic/tests/test_kernel_backends.py index eae9d2adc..b7e56cb47 100644 --- a/src/struphy/pic/tests/test_kernel_backends.py +++ b/src/struphy/pic/tests/test_kernel_backends.py @@ -194,8 +194,8 @@ def test_catalog_signatures(): for catalog in CATALOGS.values(): catalog.check_signatures() assert {package.split(".", 1)[1]: len(catalog) for package, catalog in CATALOGS.items()} == { - "pic.pushing.kernels": 45, - "pic.accumulation.kernels": 27, + "pic.pushing.kernels": 46, + "pic.accumulation.kernels": 29, "pic.diagnostics.kernels": 10, "pic.sph.kernels": 4, "pic.kernels": 4, diff --git a/src/struphy/propagators/__init__.py b/src/struphy/propagators/__init__.py index 770e51e8b..c9ba267c1 100644 --- a/src/struphy/propagators/__init__.py +++ b/src/struphy/propagators/__init__.py @@ -11,6 +11,7 @@ "CurrentCoupling5DGradB": "struphy.propagators.current_coupling_5d_gradb", "CurrentCoupling6DCurrent": "struphy.propagators.current_coupling_6d_current", "CurrentCoupling6DDensity": "struphy.propagators.current_coupling_6d_density", + "DriftKineticMaxwellCoupling": "struphy.propagators.drift_kinetic_maxwell_coupling", "EfieldWeightsCoupling": "struphy.propagators.efield_weights_coupling", "FaradayExtended": "struphy.propagators.faraday_extended", "GyrokineticPoissonSolve": "struphy.propagators.gyrokinetic_poisson_solve", @@ -62,6 +63,7 @@ from struphy.propagators.current_coupling_5d_gradb import CurrentCoupling5DGradB from struphy.propagators.current_coupling_6d_current import CurrentCoupling6DCurrent from struphy.propagators.current_coupling_6d_density import CurrentCoupling6DDensity + from struphy.propagators.drift_kinetic_maxwell_coupling import DriftKineticMaxwellCoupling from struphy.propagators.efield_weights_coupling import EfieldWeightsCoupling from struphy.propagators.faraday_extended import FaradayExtended from struphy.propagators.gyrokinetic_poisson_solve import GyrokineticPoissonSolve @@ -114,6 +116,7 @@ "CurrentCoupling5DGradB", "CurrentCoupling6DCurrent", "CurrentCoupling6DDensity", + "DriftKineticMaxwellCoupling", "EfieldWeightsCoupling", "FaradayExtended", "GyrokineticPoissonSolve", diff --git a/src/struphy/propagators/drift_kinetic_maxwell_coupling.py b/src/struphy/propagators/drift_kinetic_maxwell_coupling.py new file mode 100644 index 000000000..406238372 --- /dev/null +++ b/src/struphy/propagators/drift_kinetic_maxwell_coupling.py @@ -0,0 +1,697 @@ +"Particle and FEEC variables are updated." + +import logging +from dataclasses import dataclass +from typing import Literal + +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.dk_maxwell import dk_maxwell +from struphy.pic.accumulation.kernels.dk_maxwell_mf import dk_maxwell_mf +from struphy.pic.accumulation.kernels.marker_operator_asym_dot import marker_operator_asym_dot +from struphy.pic.accumulation.kernels.marker_operator_rank1_dot import marker_operator_rank1_dot +from struphy.pic.accumulation.particles_to_grid import Accumulator, AccumulatorVector, MarkerOperator +from struphy.pic.pushing.kernels.push_vpar_with_efield import push_vpar_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 DriftKineticMaxwellCoupling(Propagator): + r"""PIC-FEEC discretization of the following equations: + find :math:`\mathbf E \in H(\textnormal{curl})`, :math:`\tilde{\mathbf B} \in H(\textnormal{div})` and the + guiding-center distribution :math:`f(\mathbf X, v_\parallel, \mu)` (positions fixed in this step) such that + + .. math:: + + \lambda^2 \int_{\Omega} \frac{\partial \mathbf E}{\partial t} \cdot \mathbf F \,\textrm d \mathbf x + + \chi \int_{\Omega} \frac{n_0}{|B_0|^2} \frac{\partial \mathbf E_\perp}{\partial t} \cdot \mathbf F_\perp + \,\textrm d \mathbf x + &= c_1 \int_{\Omega} \tilde{\mathbf B} \cdot \nabla \times \mathbf F \,\textrm d \mathbf x + - c_2 \int_{\Omega} \iint f^\textnormal{vol} \left( v_\parallel \frac{\mathbf B^*}{B^*_\parallel} + + \frac{\mathbf E \times \mathbf b_0}{B^*_\parallel} \right) \cdot \mathbf F + \,\textrm d v_\parallel \textrm d \mu \,\textrm d \mathbf x + \qquad \forall \, \mathbf F \in H(\textnormal{curl})\,, + \\[2mm] + \frac{\partial \tilde{\mathbf B}}{\partial t} + \nabla\times\mathbf E &= 0\,, + \\[2mm] + \frac{\partial f}{\partial t} + c_3\, \frac{\mathbf B^*}{B^*_\parallel} \cdot \mathbf E\, + \frac{\partial f}{\partial v_\parallel} &= 0\,, + + with :math:`\mathbf E_\perp = \mathbf E - (\mathbf b_0 \cdot \mathbf E) \mathbf b_0`, the phase-space + density :math:`f^\textnormal{vol} = B^*_\parallel f` (the marker weights) and, as in + :class:`~struphy.propagators.push_guiding_center_parallel.PushGuidingCenterParallel`, + + .. math:: + + \mathbf B^* = \mathbf B_0 + \varepsilon v_\parallel \nabla \times \mathbf b_0\,,\qquad + B^*_\parallel = \mathbf b_0 \cdot (\mathbf B_0 + \tilde{\mathbf B}) + \varepsilon v_\parallel + (\nabla \times \mathbf b_0) \cdot \mathbf b_0\,. + + The current is the guiding-center current :math:`\iint f^\textnormal{vol} \dot{\mathbf X}` without the + magnetic drifts (which move the positions): parallel streaming, including the curvature drift in + :math:`\mathbf B^*`, and (option ``exb_current``) the :math:`E \times B` drift. The coefficients are + :math:`\lambda \geq 0`, :math:`\chi \geq 0` with :math:`\lambda^2 + \chi > 0`, :math:`c_1 > 0` and + :math:`c_2 c_3 > 0`. The defaults are the coefficients of the species' :ref:`equation parameters `, + :math:`\lambda = \varepsilon \lambda_\textrm{s}`, :math:`\chi = (\varepsilon \kappa_\textrm{s})^2`, + :math:`c_1 = 1`, :math:`c_2 = \kappa_\textrm{s}^2 \varepsilon`, :math:`c_3 = 1/\varepsilon`, i.e. for the thermal + velocity scale of :class:`~struphy.models.drift_kinetic_maxwell.DriftKineticMaxwell` + :math:`\lambda = \hat v_\textrm{th}/c`, :math:`\chi = \hat v_\textrm{th}^2/\hat v_A^2`. The :math:`\chi` term is the + (linearized, long-wavelength) polarization current of the guiding centers, which the drift-kinetic ordering + removes from :math:`\dot{\mathbf X}`: without it, the perpendicular electrostatic field is undetermined for + :math:`\lambda = 0` (flute modes), and with only the displacement current the shear Alfvén wave travels at the + speed of light. + + 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:`\tilde{\mathbf B} \in V^2_h`, and :math:`N` markers with logical + positions :math:`\boldsymbol \eta_p` (fixed), parallel velocities :math:`v_{\parallel,p}` and weights :math:`w_p`, + the semi-discrete system reads + + .. math:: + + \mathbb M_E \dot{\mathbf e} = c_1\, \mathbb C^\top \mathbb M_2 \mathbf b + - c_2 \left( \mathbf j(\mathbf V) + \mathbb S\, \mathbf e \right) \,, + \qquad + \dot{\mathbf b} = - \mathbb C \mathbf e\,, + \qquad + \dot v_{\parallel,p} = c_3\, \mathbf a_p \cdot \mathbf \Lambda^1_p \mathbf e\,, + + with :math:`\mathbb M_E = \lambda^2\, \mathbb M_1 + \chi\, \mathbb M_1^\perp` (:math:`\mathbb M_1^\perp` = + :attr:`~struphy.feec.mass.WeightedMassOperators.M1gyro`), :math:`\mathbf a_p = + \hat{\mathbf B}^{*2} / \hat B^{*3}_\parallel (\boldsymbol \eta_p)`, :math:`\mathbf \Lambda^1_p` the basis of + :math:`V^1_h` at :math:`\boldsymbol \eta_p`, and + + .. math:: + + \mathbf j(\mathbf V) = \sum_p w_p v_{\parallel,p}\, (\mathbf \Lambda^1_p)^\top \mathbf a_p\,, + \qquad + \mathbb A_\parallel = \sum_p w_p\, (\mathbf \Lambda^1_p)^\top \mathbf a_p \mathbf a_p^\top \mathbf \Lambda^1_p\,, + \qquad + \mathbf f^\top \mathbb S\, \mathbf e = \sum_p \frac{w_p}{\hat B^{*3}_\parallel(\boldsymbol \eta_p)}\, + (\mathbf \Lambda^1_p \mathbf f) \cdot \big( (\mathbf \Lambda^1_p \mathbf e) \times \hat{\mathbf b}^1_0 + (\boldsymbol \eta_p) \big)\,. + + :math:`\mathbb A_\parallel` is symmetric positive semi-definite (the parallel response), :math:`\mathbb S` is + antisymmetric (the :math:`E \times B` current does no work). The system conserves the energy + + .. math:: + + H = \frac{1}{2}\, \mathbf e^\top \mathbb M_E \mathbf e + \frac{c_1}{2}\, \mathbf b^\top \mathbb M_2 \mathbf b + + \frac{c_2}{c_3} \sum_p w_p \left( \frac{v_{\parallel,p}^2}{2} + \mu_p |B_0|(\boldsymbol \eta_p) \right) + \qquad (:meth:`~struphy.propagators.drift_kinetic_maxwell_coupling.DriftKineticMaxwellCoupling.energy`)\,. + + :ref:`time_discret`: Crank-Nicolson (implicit mid-point), with :math:`\mathbf a_p` and :math:`\mathbb S` + evaluated at the beginning of the step, which conserves :math:`H` exactly (up to the solver tolerance; for + constant weights). Eliminating :math:`\mathbf b^\textnormal{mid}` and :math:`\mathbf V^\textnormal{mid}` gives the + Schur system for the mid-point field + + .. math:: + + \left( \mathbb M_E + \frac{\Delta t^2}{4} \left( c_1\, \mathbb C^\top \mathbb M_2 \mathbb C + + c_2 c_3\, \mathbb A_\parallel \right) + \frac{\Delta t}{2}\, c_2\, \mathbb S \right) \mathbf e^\textnormal{mid} + = \mathbb M_E \mathbf e^n + \frac{\Delta t}{2} \left( c_1\, \mathbb C^\top \mathbb M_2 \mathbf b^n + - c_2\, \mathbf j(\mathbf V^n) \right)\,, + + accumulated by :func:`~struphy.pic.accumulation.kernels.dk_maxwell.dk_maxwell_kernels.dk_maxwell` (or + matrix-free) and solved with ``solver`` and ``precond`` (non-symmetric with the :math:`E \times B` current). + Then + + .. math:: + + \mathbf b^{n+1} = \mathbf b^n - \Delta t\, \mathbb C \mathbf e^\textnormal{mid}\,, + \qquad + v_{\parallel,p}^{n+1} = v_{\parallel,p}^n + \Delta t\, c_3\, \mathbf a_p \cdot \mathbf \Lambda^1_p + \mathbf e^\textnormal{mid}\,, + \qquad + \mathbf e^{n+1} = 2 \mathbf e^\textnormal{mid} - \mathbf e^n\,, + + with the velocity update by + :func:`~struphy.pic.pushing.kernels.push_vpar_with_efield.push_vpar_with_efield_kernels.push_vpar_with_efield`. + + **Quasi-neutral limit** :math:`\lambda = 0`. The parallel part of :math:`\mathbf e` is then determined by + :math:`\mathbb A_\parallel` alone (which needs positive weights, i.e. full-f, and markers everywhere), the + perpendicular part by :math:`\chi\, \mathbb M_1^\perp` and the curl-curl operator. The perpendicular field + still carries the energy :math:`\frac 12 \mathbf e^\top \mathbb M_E \mathbf e`, so the Crank-Nicolson update is + kept (storing :math:`\mathbf e^\textnormal{mid}` would dissipate it). The parallel part of the stored + :math:`\mathbf e^{n+1}` is not determined by the scheme and may oscillate from step to step; it does not feed + back, since the steps see :math:`\mathbf e^n` only through :math:`\mathbb M_E \mathbf e^n` and the + :math:`E \times B` push only through :math:`\mathbf E_\perp`. With ``initial_projection`` + (default), the setup replaces :math:`(\mathbf b^0, \mathbf V^0)` by the mid-point state of one step, which + removes the undamped period-2 oscillation of the constraint residual. + + **Preconditioner.** With :math:`\mathbb A_\parallel \approx \bar n\, \mathbb M_1^\parallel` + (:attr:`~struphy.feec.mass.WeightedMassOperators.M1para`), the symmetric part of the Schur operator is + approximately :math:`\mathbb C^\top \left(\frac{\Delta t^2 c_1}{4} \mathbb M_2\right) \mathbb C + \mathbb M_E + + \frac{\Delta t^2 c_2 c_3 \bar n}{4} \mathbb M_1^\parallel`, preconditioned (default) by + :class:`~struphy.feec.preconditioner.StiffnessPreconditioner`, built for the current time step + (:meth:`~struphy.propagators.drift_kinetic_maxwell_coupling.DriftKineticMaxwellCoupling.prepare_time_step`). + """ + + class Variables: + """Container for variables advanced by :class:`~struphy.propagators.drift_kinetic_maxwell_coupling.DriftKineticMaxwellCoupling`. + + Attributes + ---------- + e : FEECVariable + Electric-field variable in ``"Hcurl"`` space. + b : FEECVariable + Magnetic-field perturbation in ``"Hdiv"`` space. + ions : PICVariable + Guiding-center particle variable in ``"Particles5D"`` 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 == "Particles5D" + self._ions = new + + def __init__(self): + self.variables = self.Variables() + + @dataclass(repr=False) + class Options(OptionsBase): + r"""Configuration options for :class:`~struphy.propagators.drift_kinetic_maxwell_coupling.DriftKineticMaxwellCoupling`. + + Parameters + ---------- + lam : float, default=None + The displacement current has the coefficient :math:`\lambda^2` (``lam`` :math:`= \lambda \geq 0`); + ``lam=0.0`` is the quasi-neutral limit (needs ``chi > 0`` and full-f weights, no control variate). If + ``None``, :math:`\varepsilon \lambda_\textrm{s}` of the species (``lam`` and ``epsilon`` of + its equation parameters). + + chi : float, default=None + Coefficient :math:`\chi \geq 0` of the polarization current + :math:`\frac{n_0}{|B_0|^2} \partial_t \mathbf E_\perp`. If ``None``, + :math:`(\varepsilon \kappa_\textrm{s})^2` of the species (the polarization of the + guiding centers themselves); a larger value mimics the polarization of a heavier species. ``lam**2 + chi`` + must be positive. The density :math:`n_0` is the one of the MHD equilibrium + (:attr:`~struphy.feec.mass.WeightedMassOperators.M1gyro`), not of the kinetic background: choose them + consistently. + + c1 : float, default=1.0 + Coefficient :math:`c_1 > 0` of :math:`\nabla \times \tilde{\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 parallel electric force (:math:`c_2 c_3 > 0`). If ``None``, :math:`1/\varepsilon` of + the species. + + exb_current : bool, default=True + Include the :math:`E \times B` current (antisymmetric operator :math:`\mathbb S`). The Schur system is then + non-symmetric and needs a general solver. Switch it off if the :math:`E \times B` current is compensated + (e.g. by other species). + + solver : OptsSolver, default="pbicgstab" + Iterative solver of the Schur system, passed to ``feectools.linalg.solvers.inverse``. The symmetric + solvers ``"pcg"`` and ``"cg"`` need ``exb_current=False``. + + precond : LiteralOptions.OptsMaxwellPrecond, default="StiffnessPreconditioner" + Preconditioner of the Schur system (only for ``"pcg"`` and ``"pbicgstab"``; set to ``None`` otherwise). + + - ``"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`. + - ``None``: no preconditioner. + + 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_\parallel \mathbf 1 / \mathbf 1^\top \mathbb M_1^\parallel \mathbf 1` at + the setup (a density average). + + matrix_free : bool, default=False + Apply the particle matrices :math:`\mathbb A_\parallel` and :math:`\mathbb S` matrix-free, with their + fillings stored per marker (:class:`~struphy.pic.accumulation.particles_to_grid.MarkerOperator`), instead + of assembling them. Needs 7 auxiliary marker columns. + + 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 the constraint of the + quasi-neutral limit exactly. Not applied to a restored (restart) state. + """ + + # specific literals + OptsSolver = Literal["pbicgstab", "bicgstab", "gmres", "pcg", "cg"] + + lam: float = None + chi: float = None + c1: float = 1.0 + c2: float = None + c3: float = None + exb_current: bool = True + solver: OptsSolver = "pbicgstab" + 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, self.OptsSolver) + 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}." + if self.chi is not None: + assert self.chi >= 0.0, f"chi must be >= 0, got {self.chi}." + if self.lam is not None and self.chi is not None: + assert self.lam**2 + self.chi > 0.0, ( + "lam**2 + chi must be > 0: without displacement and polarization current, the perpendicular " + "electric field is undetermined." + ) + 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.exb_current, bool) + assert isinstance(self.matrix_free, bool) + assert isinstance(self.initial_projection, bool) + if self.exb_current: + assert self.solver not in ("pcg", "cg"), ( + f'solver="{self.solver}" needs a symmetric system: use a general solver or exb_current=False.' + ) + if self.solver not in ("pcg", "pbicgstab") and self.precond is not None: + logger.warning(f'solver="{self.solver}" 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 + self._epsilon = epsilon + # lam^2 = (epsilon lam_s)^2, chi = (epsilon kappa_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._chi = (epsilon * equation_params.kappa) ** 2 if self.options.chi is None else self.options.chi + 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 + self._exb = self.options.exb_current + assert self._c2 * self._c3 > 0.0, f"c2 * c3 must be > 0, got c2={self._c2}, c3={self._c3}." + assert self._lam2 + self._chi > 0.0, f"lam**2 + chi must be > 0, got lam**2={self._lam2}, chi={self._chi}." + + assert not self.derham.polar_splines, "DriftKineticMaxwellCoupling 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( + "DriftKineticMaxwellCoupling with control variate: the energy is not conserved exactly " + "(the weights change), and the response of the background f0 (parallel response and E x B current) " + "is missing in the particle matrices." + ) + + self._info = self.options.solver_params.info + + # magnetic background: B* (2-form, equilibrium) and B*_parallel with the perturbation b (frozen in a step) + self._b2 = self.projected_equil.b2 + self._curl_unit_b2 = self.projected_equil.curl_unit_b2 + self._unit_b1 = self.projected_equil.unit_b1 + self._absB0 = self.projected_equil.absB0 + self._curl_unit_b_dot_b0 = self.projected_equil.curl_unit_b_dot_b0 + + self._E0T = self.derham.extraction_ops["0"].transpose() + self._PB = getattr(self.basis_ops, "PB") + self._B_dot_b = self.derham.V0.zeros() + self._B_dot_b_pol = self._PB.codomain.zeros() + self._update_b_dot_b() + + args_fields = ( + self._b2[0]._data, + self._b2[1]._data, + self._b2[2]._data, + self._curl_unit_b2[0]._data, + self._curl_unit_b2[1]._data, + self._curl_unit_b2[2]._data, + ) + args_unit_b = ( + self._unit_b1[0]._data, + self._unit_b1[1]._data, + self._unit_b1[2]._data, + ) + args_b_par = ( + self._B_dot_b._data, + self._curl_unit_b_dot_b0._data, + ) + + # accumulation of A_par and S (assembled, or matrix-free) and the vector j + assert xp.get_backend() != "cupy", "DriftKineticMaxwellCoupling has no CUDA kernels yet." + self._matrix_free = self.options.matrix_free + if self._matrix_free: + # the vector, and the fillings of the matrices per marker (in the first free columns of the marker array) + column = particles.args_markers.first_free_idx + assert particles.n_cols_aux >= 7, "matrix_free needs 7 free (auxiliary) marker columns." + self._args_accum_pre = (epsilon,) + self._args_accum_post = args_fields + args_unit_b + args_b_par + (column,) + self._accum = AccumulatorVector( + particles, + "Hcurl", + dk_maxwell_mf, + self.mass_ops, + self.domain.args_domain, + ) + self._A_par = MarkerOperator( + particles, + "Hcurl", + marker_operator_rank1_dot, + self.mass_ops, + self.domain.args_domain, + (1, column), + symmetry="symm", + ) + self._S = MarkerOperator( + particles, + "Hcurl", + marker_operator_asym_dot, + self.mass_ops, + self.domain.args_domain, + (1, column + 4), + symmetry="asym", + ) + else: + # one matrix A_par + skew * S, skew = 2 / (dt * c3) (set in prepare_time_step) + self._args_accum_pre = (epsilon,) + self._args_accum_post = args_fields + args_unit_b + args_b_par + self._accum = Accumulator( + particles, + "Hcurl", + dk_maxwell, + self.mass_ops, + self.domain.args_domain, + add_vector=True, + symmetry=None, + ) + self._A_lam = self._accum.operators[0] + self._skew = 0.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 + if self._lam2 > 0.0 and self._chi > 0.0: + self._ME = self._lam2 * self._M1 + self._chi * self.mass_ops.M1gyro + elif self._lam2 > 0.0: + self._ME = self._lam2 * self._M1 + else: + self._ME = self._chi * self.mass_ops.M1gyro + + # reference density for the stiffness preconditioner: density average of the parallel particle matrix + if self.options.n_ref is None: + self._accumulate() + ones = self.derham.V1.zeros() + for block in ones.blocks: + block._data[:] = 1.0 + A_par = self._A_par if self._matrix_free else self._A_lam + ones_A = ones.inner(A_par.dot(ones)) + ones_M = ones.inner(self.mass_ops.M1para.dot(ones)) + if ones_A > 0.0 and ones_M > 0.0: + self._n_ref = ones_A / ones_M + else: + self._n_ref = 1.0 + logger.warning( + "DriftKineticMaxwellCoupling: 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 * a_p . Lambda^1 e_mid + args_kernel = ( + ( + self.derham.args_derham, + epsilon, + self._e_mid.blocks[0]._data, + self._e_mid.blocks[1]._data, + self._e_mid.blocks[2]._data, + ) + + args_fields + + args_b_par + + (self._c3,) + ) + + self._pusher = Pusher( + particles, + push_vpar_with_efield, + args_kernel, + self.domain.args_domain, + alpha_in_kernel=1.0, + pushes_eta=False, + ) + + def _update_b_dot_b(self): + """Parallel magnetic field b0 . (B0 + b) in the tensor-product space V0 (used in B*_parallel).""" + B_dot_b_pol = self._PB.dot(self.variables.b.spline.vector, out=self._B_dot_b_pol) + B_dot_b = self._E0T.dot(B_dot_b_pol, out=self._B_dot_b) + B_dot_b += self._absB0 + B_dot_b.update_ghost_regions() + + def _accumulate(self): + """Accumulate the particle matrices (or their fillings per marker) and the vector j at the current markers.""" + if self._matrix_free: + self._accum(*self._args_accum_pre, *self._args_accum_post) + else: + self._accum(*self._args_accum_pre, self._skew, *self._args_accum_post) + + def prepare_time_step(self, dt: float): + """Schur operator, its preconditioner and solver for the time step ``dt``: built in the setup, rebuilt when dt + changes. The particle matrices enter by reference (assembled, or matrix-free reading the fillings stored in + the markers), so they are updated by the accumulation in every step. For ``lam=0`` (and + ``initial_projection``), the first call also projects the state onto the constraint.""" + if dt == self._schur_dt: + return + + schur = self._ME + (dt**2 * self._c1 / 4) * self._curl_curl + if self._matrix_free: + schur = schur + (dt**2 * self._c2 * self._c3 / 4) * self._A_par + if self._exb: + schur = schur + (dt * self._c2 / 2) * self._S + else: + # A_lam = A_par + skew S with skew = 2 / (dt c3): dt^2 c2 c3 / 4 A_lam = dt^2 c2 c3 / 4 A_par + dt c2 / 2 S + self._skew = 2.0 / (dt * self._c3) if self._exb else 0.0 + schur = schur + (dt**2 * self._c2 * self._c3 / 4) * self._A_lam + + if self._precond == "StiffnessPreconditioner": + stab = self._ME + (dt**2 * self._c2 * self._c3 * self._n_ref / 4) * self.mass_ops.M1para + pc = preconditioner.StiffnessPreconditioner( + self.mass_ops, + "curl", + sigma=1.0, + mid=(dt**2 * self._c1 / 4) * self._M2, + stab=stab, + ) + 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{1}{2} \mathbf e^\top \mathbb M_E \mathbf e` (displacement and polarization), + :math:`\frac{c_1}{2} \mathbf b^\top \mathbb M_2 \mathbf b` and + :math:`\frac{c_2}{c_3} \sum_p w_p (v_{\parallel,p}^2/2 + \mu_p |B_0|(\boldsymbol \eta_p))` 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 + + particles.save_magnetic_background_energy() + valid = particles.valid_mks + markers = particles.markers + kinetic = xp.array( + [ + xp.sum( + markers[5, valid] * (0.5 * markers[3, valid] ** 2 + markers[particles.first_diagnostics_idx, valid]) + ) + ] + ) + 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 = 0.5 * e.inner(self._ME.dot(e, out=self._e_tmp)) + 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{1}{2} \mathbf e^\top \mathbb M_E \mathbf e + \frac{c_1}{2} \mathbf b^\top + \mathbb M_2 \mathbf b + \frac{c_2}{c_3} \sum_p w_p (v_{\parallel,p}^2/2 + \mu_p |B_0|(\boldsymbol \eta_p))` + 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 + + # B*_parallel with b^n, then A_par, S (or their fillings per marker) and j + self._update_b_dot_b() + self._accumulate() + + # rhs = M_E 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 + rhs += self._ME.dot(en, out=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 the constraint of the + quasi-neutral limit exactly. The electric field is kept.""" + 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 a_p . 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("DriftKineticMaxwellCoupling: projected the initial state onto the constraint (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 (M_E = lam^2 M1 + chi M1gyro is nonzero, so e carries energy also for lam = 0) + e_new = e_mid.copy(out=self._e_new) + e_new *= 2.0 + e_new -= en + + # update parallel 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 DriftKineticMaxwellCoupling: {info['success']}") + logger.info(f"Iterations for DriftKineticMaxwellCoupling: {info['niter']}") + logger.info(f"Maxdiff e for DriftKineticMaxwellCoupling: {diffs['e']}") + logger.info(f"Maxdiff b for DriftKineticMaxwellCoupling: {diffs['b']}") + logger.info("") diff --git a/src/struphy/propagators/push_guiding_center_bx_estar.py b/src/struphy/propagators/push_guiding_center_bx_estar.py index f6fc1da6f..c1a08dde9 100644 --- a/src/struphy/propagators/push_guiding_center_bx_estar.py +++ b/src/struphy/propagators/push_guiding_center_bx_estar.py @@ -1,7 +1,7 @@ "Only particle variables are updated." import logging -from dataclasses import dataclass +from dataclasses import dataclass, replace from typing import Literal from line_profiler import profile @@ -42,7 +42,7 @@ class PushGuidingCenterBxEstar(Propagator): \mathbf E^* = -\nabla \phi - \varepsilon \mu_p \nabla |\mathbf B|\,,\qquad \mathbf B^* = \mathbf B + \varepsilon v_\parallel \nabla \times \mathbf{b}_0\,,\qquad B^*_\parallel = \mathbf B^* \cdot \mathbf{b}_0\,, - Here, :math:`\mathbf B = \mathbf{B}_0 + \tilde{\mathbf B}` can be the full magnetic field (equilibrium + perturbation). The electric potential :math:`\phi` and/or the magnetic perturbation :math:`\tilde{\mathbf B}` can be ignored by passing ``None``. + Here, :math:`\mathbf B = \mathbf{B}_0 + \tilde{\mathbf B}` can be the full magnetic field (equilibrium + perturbation). The electric potential :math:`\phi` and/or the magnetic perturbation :math:`\tilde{\mathbf B}` can be ignored by passing ``None``. Instead of the potential, an electric field :math:`\mathbf E \in H(\textnormal{curl})` can be passed (``e_field``, explicit algorithm); then :math:`\mathbf E^* = \mathbf E - \varepsilon \mu_p \nabla |\mathbf B|` (gauge-invariant electromagnetic formulation). In logical space this is given by :math:`\mathbf X = F(\boldsymbol \eta)`: .. math:: @@ -83,8 +83,9 @@ def __init__( self, phi: FEECVariable = None, b_tilde: FEECVariable = None, + e_field: FEECVariable = None, ): - """ + r""" Parameters ---------- phi : FEECVariable, default=None @@ -93,7 +94,16 @@ def __init__( b_tilde : FEECVariable, default=None Magnetic perturbation in ``"Hcurl"`` space contributing to :math:`B^*`. If ``None``, the perturbation is ignored. + e_field : FEECVariable, default=None + Electric field in ``"Hcurl"`` space, used in :math:`E^* = E - \varepsilon \mu \nabla |B|` instead of + :math:`-\nabla \phi` (gauge-invariant electromagnetic formulation). Cannot be combined with ``phi``. Only the + explicit algorithm is consistent with a general :math:`E` (the discrete gradients conserve a Hamiltonian + with a potential); a discrete-gradient ``algo`` is replaced by ``"explicit"``, with a warning. """ + assert phi is None or e_field is None, "Pass either the potential phi or the electric field e_field, not both." + if e_field is not None: + assert isinstance(e_field, FEECVariable) and e_field.space == "Hcurl", "e_field must be in Hcurl." + self.e_field = e_field self.variables = self.Variables() self.phi = phi if phi is not None else FEECVariable(space="H1") if phi is None: @@ -203,6 +213,16 @@ def allocate(self): self._phi = self.derham.V0.zeros() self._update_phi() self._evaluate_e_field = self.options.evaluate_e_field + if self.e_field is not None: + self._evaluate_e_field = True + # the discrete gradients conserve the Hamiltonian, which contains no potential for a general E: + # they would remove the work of E along the step (inconsistent scheme), hence the explicit algorithm + if "discrete_gradient" in self.options.algo: + logger.warning( + f'{self.__class__.__name__} with e_field: algo="{self.options.algo}" needs a potential, ' + 'using algo="explicit" (default Butcher tableau).' + ) + self._options = replace(self.options, algo="explicit", butcher=None) self._e_field = self.derham.V1.zeros() self._e_field_pol = self.derham.grad.codomain.zeros() @@ -509,7 +529,10 @@ def __call__(self, dt): # electric field # TODO: add out to __neg__ of StencilVector - if self._evaluate_e_field: + if self.e_field is not None: + e_field = self._E1T.dot(self.e_field.spline.vector, out=self._e_field) + e_field.update_ghost_regions() + elif self._evaluate_e_field: e_field_pol = self.derham.grad.dot(-self.phi.spline.vector, out=self._e_field_pol) e_field = self._E1T.dot(e_field_pol, out=self._e_field) e_field.update_ghost_regions() diff --git a/src/struphy/propagators/push_guiding_center_parallel.py b/src/struphy/propagators/push_guiding_center_parallel.py index 7aa53f662..8c6472961 100644 --- a/src/struphy/propagators/push_guiding_center_parallel.py +++ b/src/struphy/propagators/push_guiding_center_parallel.py @@ -1,7 +1,7 @@ "Only particle variables are updated." import logging -from dataclasses import dataclass +from dataclasses import dataclass, replace from typing import Literal from line_profiler import profile @@ -52,7 +52,7 @@ class PushGuidingCenterParallel(Propagator): where :math:`\mathbf B = \mathbf B_0 + \tilde{\mathbf B}` can be the full magnetic field (equilibrium + perturbation). The electric potential ``phi`` and/or the magnetic perturbation ``b_tilde`` - can be ignored by passing ``None``. + can be ignored by passing ``None``. Instead of the potential, an electric field :math:`\mathbf E \in H(\textnormal{curl})` can be passed (``e_field``, explicit algorithm); then :math:`\mathbf E^* = \mathbf E - \varepsilon \mu_p \nabla |\mathbf B|` (gauge-invariant electromagnetic formulation). In logical space this is given by :math:`\mathbf X = F(\boldsymbol \eta)`: .. math:: @@ -99,8 +99,9 @@ def __init__( self, phi: FEECVariable = None, b_tilde: FEECVariable = None, + e_field: FEECVariable = None, ): - """ + r""" Parameters ---------- phi : FEECVariable, default=None @@ -109,7 +110,16 @@ def __init__( b_tilde : FEECVariable, default=None Magnetic perturbation in ``"Hcurl"`` space contributing to :math:`B^*`. If ``None``, the perturbation is ignored. + e_field : FEECVariable, default=None + Electric field in ``"Hcurl"`` space, used in :math:`E^* = E - \varepsilon \mu \nabla |B|` instead of + :math:`-\nabla \phi` (gauge-invariant electromagnetic formulation). Cannot be combined with ``phi``. Only the + explicit algorithm is consistent with a general :math:`E` (the discrete gradients conserve a Hamiltonian + with a potential); a discrete-gradient ``algo`` is replaced by ``"explicit"``, with a warning. """ + assert phi is None or e_field is None, "Pass either the potential phi or the electric field e_field, not both." + if e_field is not None: + assert isinstance(e_field, FEECVariable) and e_field.space == "Hcurl", "e_field must be in Hcurl." + self.e_field = e_field self.variables = self.Variables() self.phi = phi if phi is not None else FEECVariable(space="H1") if phi is None: @@ -220,6 +230,16 @@ def allocate(self): self._phi = self.derham.V0.zeros() self._update_phi() self._evaluate_e_field = self.options.evaluate_e_field + if self.e_field is not None: + self._evaluate_e_field = True + # the discrete gradients conserve the Hamiltonian, which contains no potential for a general E: + # they would remove the work of E along the step (inconsistent scheme), hence the explicit algorithm + if "discrete_gradient" in self.options.algo: + logger.warning( + f'{self.__class__.__name__} with e_field: algo="{self.options.algo}" needs a potential, ' + 'using algo="explicit" (default Butcher tableau).' + ) + self._options = replace(self.options, algo="explicit", butcher=None) self._e_field = self.derham.V1.zeros() self._e_field_pol = self.derham.grad.codomain.zeros() @@ -535,7 +555,10 @@ def __call__(self, dt): # electric field # TODO: add out to __neg__ of StencilVector - if self._evaluate_e_field: + if self.e_field is not None: + e_field = self._E1T.dot(self.e_field.spline.vector, out=self._e_field) + e_field.update_ghost_regions() + elif self._evaluate_e_field: e_field_pol = self.derham.grad.dot(-self.phi.spline.vector, out=self._e_field_pol) e_field = self._E1T.dot(e_field_pol, out=self._e_field) e_field.update_ghost_regions() diff --git a/src/struphy/propagators/tests/test_drift_kinetic_maxwell_coupling.py b/src/struphy/propagators/tests/test_drift_kinetic_maxwell_coupling.py new file mode 100644 index 000000000..1f0db84dc --- /dev/null +++ b/src/struphy/propagators/tests/test_drift_kinetic_maxwell_coupling.py @@ -0,0 +1,305 @@ +import logging + +import cunumpy as xp +import pytest + +from struphy import set_logging_level + +set_logging_level(logging.WARNING) + + +def _model(control_variate=False, domain=None, **options): + """Allocated DriftKineticMaxwell model (slab with a tilted homogeneous B0, full-f guiding centers with epsilon=0.5, + lam=2, kappa=2, a parallel flow and a B perturbation) with the given options of its coupling propagator.""" + from struphy import ( + DerhamOptions, + LoadingParameters, + Simulation, + WeightsParameters, + domains, + equils, + grids, + maxwellians, + models, + perturbations, + ) + from struphy.linear_algebra.solver import SolverParameters + + model = models.DriftKineticMaxwell(epsilon=0.5, lam=2.0, kappa=2.0) + model.kinetic_ions.set_markers( + loading_params=LoadingParameters(ppc=20, seed=1234), + weights_params=WeightsParameters(control_variate=control_variate), + ) + model.kinetic_ions.var.add_background(maxwellians.GyroMaxwellian2D(n=(1.0, None), B0=1.0)) + model.kinetic_ions.var.add_initial_condition( + maxwellians.GyroMaxwellian2D( + n=(1.0, None), u_para=(0.0, perturbations.ModesCos(ls=(1,), ns=(1,), amps=(1e-1,))), B0=1.0 + ) + ) + model.em_fields.b_field.add_perturbation( + perturbations.ModesSin(ls=(1,), ns=(1,), given_in_basis="2", comp=1, amps=(1e-1,)) + ) + + prop = model.propagators.coupling + prop.options = prop.Options(solver_params=SolverParameters(tol=1e-13, maxiter=3000), **options) + + Simulation( + model=model, + domain=domains.Cuboid(r1=4.0, r2=2.0, r3=4.0) if domain is None else domain, + equil=equils.HomogenSlab(B0x=0.3, B0z=1.0), + grid=grids.TensorProductGrid(num_elements=(6, 2, 6)), + derham_opts=DerhamOptions(degree=(2, 1, 2)), + ).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.markers[3, model.kinetic_ions.var.particles.valid_mks].copy() + return e, b, v + + +def _norm(x): + return x.inner(x) ** 0.5 + + +def _random_v1(derham, seed): + """A random V1 coefficient vector (different on each rank, which does not matter for the tests).""" + rng = xp.random.default_rng(seed) + x = derham.V1.zeros() + for block in x.blocks: + block._data[:] = rng.standard_normal(block._data.shape) + x.update_ghost_regions() + return x + + +@pytest.mark.parametrize("matrix_free", [False, True]) +@pytest.mark.parametrize("exb_current", [True, False]) +@pytest.mark.parametrize("lam", [1.0, 0.0]) +def test_energy_conservation(lam, exb_current, matrix_free): + """The coupling step conserves its energy H exactly (up to the solver tolerance), also for a large time step and + in the quasi-neutral limit lam=0 (without the initial projection, which changes H once).""" + solver = {} if exb_current else {"solver": "pcg"} + model = _model(lam=lam, exb_current=exb_current, matrix_free=matrix_free, initial_projection=False, **solver) + 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(0.5) + assert prop._solver._info["success"] + # H is dominated by the thermal energy: compare its change with the energy exchanged with the B field + exchanged = abs(0.5 * b.inner(prop._M2.dot(b)) - energy_b_0) + assert exchanged > 0.1 * energy_b_0 + assert abs(prop.energy() - energy_0) < 1e-10 * exchanged + + +def test_energy_conservation_mapped(): + """Also in a non-trivial mapping (metric in the parallel response and in the E x B current).""" + from struphy import domains + + model = _model(domain=domains.Colella(Lx=4.0, Ly=2.0, alpha=0.05, Lz=4.0), initial_projection=False) + prop = model.propagators.coupling + energy_0 = prop.energy() + for _ in range(3): + prop(0.5) + assert prop._solver._info["success"] + assert abs(prop.energy() - energy_0) < 1e-11 * energy_0 + + +@pytest.mark.parametrize("lam", [1.0, 0.0]) +def test_matrix_free_agrees(lam): + """The matrix-free particle matrices give the same step as the assembled ones (with the E x B current).""" + states = [] + for matrix_free in (False, True): + model = _model(lam=lam, matrix_free=matrix_free, initial_projection=False) + model.propagators.coupling(0.5) + # the two simulations have different Derham objects: compare the arrays + e, b, v = _state(model) + states.append((e.toarray(), b.toarray(), v)) + for x0, x1 in zip(*states): + assert xp.max(xp.abs(x1 - x0)) < 1e-9 * xp.max(xp.abs(x0)) + + +def test_particle_matrices(): + """The parallel response A_par is symmetric positive semi-definite, the E x B operator S antisymmetric.""" + model = _model(matrix_free=True, initial_projection=False) + prop = model.propagators.coupling + prop._accumulate() + x = _random_v1(prop.derham, 1) + y = _random_v1(prop.derham, 2) + A, S = prop._A_par, prop._S + scale = x.inner(A.dot(x)) + y.inner(A.dot(y)) + assert x.inner(A.dot(x)) >= 0.0 + assert abs(x.inner(A.dot(y)) - y.inner(A.dot(x))) < 1e-12 * scale + scale_s = abs(x.inner(S.dot(y))) + assert scale_s > 0.0 + assert abs(x.inner(S.dot(x))) < 1e-12 * scale_s + assert abs(x.inner(S.dot(y)) + y.inner(S.dot(x))) < 1e-12 * scale_s + + +@pytest.mark.parametrize("matrix_free", [False, True]) +@pytest.mark.parametrize("lam", [1.0, 0.0]) +def test_discrete_equations(lam, matrix_free): + """Without the E x B current, the new state satisfies the Crank-Nicolson equations: Ampère's law + M_E (e^{n+1} - e^n) / dt = c1 C^T M2 b_mid - c2 j(V_mid) (j linear in v_par at fixed positions and fixed + a_p = B*/B*_par, which depends on v_par only through the curvature, zero here) and Faraday's law.""" + dt = 0.5 + model = _model(lam=lam, matrix_free=matrix_free, exb_current=False, solver="pcg", initial_projection=False) + prop = model.propagators.coupling + e0, b0, _ = _state(model) + prop._accumulate() + j0 = prop._accum.vectors[0].copy() + prop(dt) + e1, b1, _ = _state(model) + prop._accumulate() + j_mid = (j0 + prop._accum.vectors[0]) * 0.5 + b_mid = (b0 + b1) * 0.5 + e_mid = (e0 + e1) * 0.5 + + curl_b = prop._curl.T.dot(prop._M2.dot(b_mid)) * prop._c1 + current = j_mid * prop._c2 + residual = curl_b - current - prop._ME.dot(e1 - e0) * (1 / 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 + + +def test_initial_projection(): + """For lam=0, without the initial projection the magnetic energy oscillates with period 2 (the constraint + residual changes sign in every step); the projection removes this.""" + + def en_b_series(projection): + model = _model(lam=0.0, initial_projection=projection) + prop = model.propagators.coupling + b = model.em_fields.b_field.spline.vector + series = [] + for _ in range(6): + prop(0.5) + series.append(0.5 * b.inner(prop._M2.dot(b))) + return xp.array(series) + + def zigzag(series): + # second differences relative to the mean: large for a period-2 oscillation + return float(xp.max(xp.abs(series[2:] - 2 * series[1:-1] + series[:-2])) / xp.mean(series)) + + assert zigzag(en_b_series(False)) > 3 * zigzag(en_b_series(True)) + + +def test_dt_change(): + """The Schur solver and its preconditioner are rebuilt only when dt changes; the energy stays conserved.""" + 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-11 * energy_0 + + +@pytest.mark.mpi_skip +def test_invalid_options(): + from struphy.propagators.drift_kinetic_maxwell_coupling import DriftKineticMaxwellCoupling + + Options = DriftKineticMaxwellCoupling.Options + with pytest.raises(AssertionError): + Options(lam=-1.0) + with pytest.raises(AssertionError): + Options(chi=-1.0) + with pytest.raises(AssertionError): + # no displacement and no polarization current: the perpendicular field is undetermined + Options(lam=0.0, chi=0.0) + with pytest.raises(AssertionError): + Options(c1=0.0) + with pytest.raises(AssertionError): + Options(c2=1.0, c3=-1.0) + with pytest.raises(AssertionError): + # the E x B current makes the system non-symmetric + Options(solver="pcg") + Options(solver="pcg", exb_current=False) + + +def test_control_variate_needs_lam(): + """lam=0 needs full-f weights (the particle matrix is indefinite with the control variate).""" + with pytest.raises(AssertionError): + _model(lam=0.0, control_variate=True) + + +def _cold_model(push_para_e_field=False): + """DriftKineticMaxwell with cold guiding centers (mu ~ 0), uniform B0 and no B perturbation (no mirror force), and + a large electric field with a parallel component. Optionally pushes with the electric field also in push_para + (forbidden in the model).""" + from struphy import DerhamOptions, LoadingParameters, Simulation, domains, equils, grids, maxwellians, models + from struphy.propagators.push_guiding_center_parallel import PushGuidingCenterParallel + + vth = 1e-6 + model = models.DriftKineticMaxwell(epsilon=0.5, lam=2.0, kappa=2.0) + model.kinetic_ions.set_markers( + loading_params=LoadingParameters(ppc=5, seed=1, moments=(0.0, 0.0, 1.0, vth), B0=1.0) + ) + model.kinetic_ions.var.add_background( + maxwellians.GyroMaxwellian2D(n=(1.0, None), vth_para=(1.0, None), vth_perp=(vth, None), B0=1.0) + ) + if push_para_e_field: + model.propagators.push_para = PushGuidingCenterParallel( + b_tilde=model.em_fields.b_field, e_field=model.em_fields.e_field + ) + model.propagators.push_para.variables.ions = model.kinetic_ions.var + Simulation( + model=model, + domain=domains.Cuboid(r1=2.0, r2=2.0, r3=2.0), + equil=equils.HomogenSlab(B0x=0.3, B0z=1.0), + grid=grids.TensorProductGrid(num_elements=(4, 4, 4)), + derham_opts=DerhamOptions(degree=(2, 2, 2)), + ).allocate() + e = model.em_fields.e_field.spline.vector + for blk in e.blocks: + blk._data[:] = 1.0 + e.update_ghost_regions() + return model + + +def _all_vpar(particles): + """Sorted parallel velocities of all valid markers on all processes (markers may change their process during a + push, so they are compared as a set).""" + from maybempi import MPI + + local = particles.markers[3, particles.valid_mks].copy() + comm = MPI.COMM_WORLD + if comm.Get_size() > 1: + local = xp.concatenate(comm.allgather(local)) + return xp.sort(local) + + +@pytest.mark.parametrize("evaluate_e_field", [False, True]) +@pytest.mark.parametrize("algo", ["discrete_gradient_1st_order", "explicit"]) +def test_parallel_force_applied_once(algo, evaluate_e_field): + """The parallel electric force b* . E / epsilon is applied only by the coupling step: push_para of the model + leaves v_par unchanged when E != 0 (cold markers, no mirror force), for every algorithm. An electric field passed + to push_para in the model is rejected.""" + model = _cold_model() + prop = model.propagators.push_para + prop.options = prop.Options(algo=algo, evaluate_e_field=evaluate_e_field) + prop.allocate() + particles = model.kinetic_ions.var.particles + v_0 = _all_vpar(particles) + prop(0.1) + assert xp.max(xp.abs(_all_vpar(particles) - v_0)) == 0.0 + + # the same push with the electric field changes v_par (by about dt / epsilon * b* . E) + model.propagators.push_para.e_field = model.em_fields.e_field + prop.options = prop.Options(algo="explicit") + prop.allocate() + prop(0.1) + assert xp.max(xp.abs(_all_vpar(particles) - v_0)) > 1e-2 + + +def test_parallel_force_in_push_para_rejected(): + with pytest.raises(AssertionError): + _cold_model(push_para_e_field=True) diff --git a/src/struphy/propagators/tests/test_push_guiding_center_e_field.py b/src/struphy/propagators/tests/test_push_guiding_center_e_field.py new file mode 100644 index 000000000..0389c476e --- /dev/null +++ b/src/struphy/propagators/tests/test_push_guiding_center_e_field.py @@ -0,0 +1,98 @@ +import logging + +import cunumpy as xp +import pytest + +from struphy import set_logging_level + +set_logging_level(logging.WARNING) + + +def _allocate(model, perturb_density): + """Same markers for all models (same loading, seed, domain and grid); a density perturbation makes the initial + potential of the electrostatic model non-trivial.""" + from struphy import DerhamOptions, LoadingParameters, Simulation, domains, equils, grids, maxwellians, perturbations + + model.kinetic_ions.set_markers(loading_params=LoadingParameters(ppc=10, seed=4321)) + model.kinetic_ions.var.add_background(maxwellians.GyroMaxwellian2D(n=(1.0, None), B0=1.0)) + pert = perturbations.ModesCos(ls=(1,), ms=(1,), amps=(0.3,)) if perturb_density else None + model.kinetic_ions.var.add_initial_condition(maxwellians.GyroMaxwellian2D(n=(1.0, pert), B0=1.0)) + Simulation( + model=model, + domain=domains.Cuboid(r1=4.0, r2=4.0, r3=2.0), + equil=equils.HomogenSlab(B0x=0.2, B0z=1.0), + grid=grids.TensorProductGrid(num_elements=(6, 6, 2)), + derham_opts=DerhamOptions(degree=(2, 2, 1)), + ).allocate() + + +def _all_markers(particles): + """Positions and parallel velocities of all valid markers on all processes, sorted lexicographically (markers + may change their process during a push).""" + from maybempi import MPI + + local = particles.markers[:4, particles.valid_mks].copy() + comm = MPI.COMM_WORLD + if comm.Get_size() > 1: + local = xp.concatenate(comm.allgather(local), axis=1) + return local[:, xp.lexsort(local[::-1])] + + +@pytest.mark.parametrize("e_algo", ["explicit", "discrete_gradient_1st_order"]) +@pytest.mark.parametrize("name", ["PushGuidingCenterBxEstar", "PushGuidingCenterParallel"]) +def test_e_field_equals_minus_grad_phi(name, e_algo): + """A guiding-center pusher with ``e_field = -grad phi`` moves the markers exactly as with the potential ``phi`` + (explicit algorithm, with the electric field evaluated from phi). With ``e_field`` a discrete-gradient algorithm + is replaced by the explicit one (it would remove the work of E along the step).""" + from struphy import models, propagators + + dt = 0.1 + attr = {"PushGuidingCenterBxEstar": "push_gc_bxe", "PushGuidingCenterParallel": "push_gc_para"}[name] + + # reference: the phi path in the electrostatic drift-kinetic model (phi from its initial Poisson solve); + # pushed before the second model is allocated (the Derham objects of propagators are class attributes) + es = models.DriftKineticElectrostaticAdiabatic(epsilon=0.5) + prop_phi = getattr(es.propagators, attr) + prop_phi.options = prop_phi.Options(algo="explicit", evaluate_e_field=True) + _allocate(es, perturb_density=True) + phi = es.em_fields.phi.spline.vector + assert phi.inner(phi) > 0.0 + e_data = [blk._data.copy() for blk in es.em_fields.phi.spline.derham.grad.dot(-phi).blocks] + p_es = es.kinetic_ions.var.particles + initial = _all_markers(p_es) + prop_phi(dt) + pushed_phi = _all_markers(p_es) + assert xp.max(xp.abs(pushed_phi - initial)) > 1e-6 + + # the e_field path, with e = -grad phi copied into the (identically distributed) Derham of the other model + em = models.DriftKineticMaxwell(epsilon=0.5) + _allocate(em, perturb_density=True) + e = em.em_fields.e_field.spline.vector + for blk, data in zip(e.blocks, e_data): + blk._data[:] = data + e.update_ghost_regions() + prop_e = getattr(propagators, name)(e_field=em.em_fields.e_field) + prop_e.variables.ions = em.kinetic_ions.var + prop_e.options = prop_e.Options(algo=e_algo) + prop_e.allocate() + assert prop_e.options.algo == "explicit" + + p_em = em.kinetic_ions.var.particles + assert xp.array_equal(initial, _all_markers(p_em)) + prop_e(dt) + pushed_e = _all_markers(p_em) + assert pushed_e.shape == pushed_phi.shape + assert xp.max(xp.abs(pushed_e - pushed_phi)) < 1e-12 + + +@pytest.mark.mpi_skip +def test_phi_and_e_field_exclusive(): + from struphy.models.variables import FEECVariable + from struphy.propagators.push_guiding_center_bx_estar import PushGuidingCenterBxEstar + from struphy.propagators.push_guiding_center_parallel import PushGuidingCenterParallel + + for cls in (PushGuidingCenterBxEstar, PushGuidingCenterParallel): + with pytest.raises(AssertionError): + cls(phi=FEECVariable(space="H1"), e_field=FEECVariable(space="Hcurl")) + with pytest.raises(AssertionError): + cls(e_field=FEECVariable(space="H1"))