Added relativistic boris push

changelog/979.feature.rst

@@ -0,0 +1 @@
+Added a realtivistic version of the explicit Boris push algorithm to `simulation.particle_integrators'.

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,75 @@ def boris_push(x, v, B, E, q, m, dt):
     v[...] = vplus + hqmdt * E
 
     x += v * dt
+
+
+def boris_push_relativistic(x, v, B, E, q, m, dt):
+    r"""
+    The explicit Boris pusher, including realtivistic 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
+    ----------
+    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.
+
+    This ends up causing the magnetic field action to be properly "centered" in
+    time, and the algorithm, being a symplectic integrator, conserves energy.
+
+    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)
+    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.si.value, t)
+    uvel_plus = uvel_minus + np.cross(uvel_prime.si.value, 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

@@ -60,7 +60,10 @@ class ParticleTracker:
     .. _`Particle Stepper Notebook`: ../notebooks/particle_stepper.ipynb
     """
 
-    integrators = {"explicit_boris": particle_integrators.boris_push}
+    integrators = {
+        "explicit_boris": particle_integrators.boris_push,
+        "explicit_boris_relativistic": particle_integrators.boris_push_relativistic,
+    }
 
     _wip_integrators = {}