Added relativistic boris push

changelog/979.feature.rst

@@ -1 +1 @@
-Added a realtivistic version of the explicit Boris push algorithm to `simulation.particle_integrators'.
+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

@@ -71,9 +71,12 @@ def boris_push(x, v, B, E, q, m, dt):
     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 realtivistic corrections.
+    The explicit Boris pusher, including relativistic corrections.
 
     Parameters
     ----------
@@ -94,48 +97,38 @@ def boris_push_relativistic(x, v, B, E, q, m, dt):
 
     Notes
     ----------
-    The Boris algorithm is the standard energy conserving algorithm for
-    particle movement in plasma physics. See [1]_ for more details, and
-    [2]_ for a nice overview.
-
-    Conceptually, the algorithm has three phases:
-
-    1. Add half the impulse from electric field.
-    2. Rotate the particle velocity about the direction of the magnetic
-       field.
-    3. Add the second half of the impulse from the electric field.
+    For the basic overview of this algorithm, see `boris_push`. This
+    version, based on [1]_, applies relativistic corrections such as
+    TODO.
 
-    This ends up causing the magnetic field action to be properly "centered" in
-    time, and the algorithm, being a symplectic integrator, conserves energy.
+    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
-    .. [2] L. Brieda, "Particle Push in Magnetic Field (Boris Method)",
-           https://www.particleincell.com/2011/vxb-rotation/
     """
-    c = constants.c.si.value
 
-    γ = 1 / np.sqrt(1 - (v / c) ** 2)
+    γ = 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)
+    γ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.si.value, t)
-    uvel_plus = uvel_minus + np.cross(uvel_prime.si.value, s)
+    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)
+    γ2 = np.sqrt(1 + (uvel_new / _c) ** 2)
 
     # Update the velocities of the particles that are being pushed
     v[...] = uvel_new / γ2

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
 
@@ -62,10 +63,12 @@ class ParticleTracker:
 
     integrators = {
         "explicit_boris": particle_integrators.boris_push,
-        "explicit_boris_relativistic": particle_integrators.boris_push_relativistic,
     }
 
-    _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)
 
@@ -110,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(