Added relativistic boris push

changelog/979.feature.rst

@@ -0,0 +1 @@
+Added a draft of a relativistic explicit Boris push algorithm to `simulation.particle_integrators'. For now, this integrator is a marked as a work in progress and will raise a warning when used.

plasmapy/simulation/particle_integrators.py

@@ -14,7 +14,7 @@
 
 def boris_push(x, v, B, E, q, m, dt):
     r"""
-    The explicit Boris pusher.
+    The explicit Boris pusher (non-relativistic)
 
     Parameters
     ----------
@@ -69,3 +69,68 @@ def boris_push(x, v, B, E, q, m, dt):
     v[...] = vplus + hqmdt * E
 
     x += v * dt
+
+
+_c = constants.c.si.value
+
+
+def boris_push_relativistic(x, v, B, E, q, m, dt):
+    r"""
+    The explicit Boris pusher, including relativistic corrections.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        particle position at full timestep, in SI (meter) units.
+    v : np.ndarray
+        particle velocity at half timestep, in SI (meter/second) units.
+    B : np.ndarray
+        magnetic field at full timestep, in SI (tesla) units.
+    E : float
+        electric field at full timestep, in SI (V/m) units.
+    q : float
+        particle charge, in SI (Coulomb) units.
+    m : float
+        particle mass, in SI (kg) units.
+    dt : float
+        timestep, in SI (second) units.
+
+    Notes
+    ----------
+    For the basic overview of this algorithm, see `boris_push`. This
+    version, based on [1]_, applies relativistic corrections such as
+    TODO.
+
+    Keep in mind that the non-relativistic version will be slightly
+    faster if you don't encounter velocities in relativistic regimes.
+
+    References
+    ----------
+    .. [1] C. K. Birdsall, A. B. Langdon, "Plasma Physics via Computer
+           Simulation", 2004, p. 58-63
+    """
+
+    γ = 1 / np.sqrt(1 - (v / _c) ** 2)
+    uvel = v * γ
+
+    uvel_minus = uvel + q * E * dt / (2 * m)
+
+    γ1 = np.sqrt(1 + (uvel_minus / _c) ** 2)
+
+    # Birdsall has a factor of c incorrect in the definiton of t?
+    # See this source: https://www.sciencedirect.com/science/article/pii/S163107211400148X
+    t = q * B * dt / (2 * γ1 * m)
+    s = 2 * t / (1 + (t * t).sum(axis=1, keepdims=True))
+
+    uvel_prime = uvel_minus + np.cross(uvel_minus, t)
+    uvel_plus = uvel_minus + np.cross(uvel_prime, s)
+    uvel_new = uvel_plus + +q * E * dt / (2 * m)
+
+    # You can show that this expression is equivalent to calculating
+    # v_new  then calculating γnew using the usual formula
+    γ2 = np.sqrt(1 + (uvel_new / _c) ** 2)
+
+    # Update the velocities of the particles that are being pushed
+    v[...] = uvel_new / γ2
+
+    x += v * dt

plasmapy/simulation/particletracker.py

@@ -6,6 +6,7 @@
 import astropy.units as u
 import numpy as np
 import scipy.interpolate as interp
+import warnings
 
 from astropy import constants
 
@@ -60,9 +61,14 @@ class ParticleTracker:
     .. _`Particle Stepper Notebook`: ../notebooks/particle_stepper.ipynb
     """
 
-    integrators = {"explicit_boris": particle_integrators.boris_push}
+    integrators = {
+        "explicit_boris": particle_integrators.boris_push,
+    }
 
-    _wip_integrators = {}
+    # Work In Progress integrators (throws a warning if used)
+    _wip_integrators = {
+        "explicit_boris_relativistic": particle_integrators.boris_push_relativistic,
+    }
 
     _all_integrators = dict(**integrators, **_wip_integrators)
 
@@ -107,6 +113,12 @@ def __init__(
         _B = np.moveaxis(self.plasma.magnetic_field.si.value, 0, -1)
         _E = np.moveaxis(self.plasma.electric_field.si.value, 0, -1)
 
+        if integrator not in self.integrators.keys():
+            warnings.warn(
+                f"The {integrator} integrator is still under development. "
+                "Use at your own risk."
+            )
+
         self.integrator = self._all_integrators[integrator]
 
         self._B_interpolator = interp.RegularGridInterpolator(