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 warn when used.

docs/notebooks/particle_stepper.ipynb

@@ -157,7 +157,9 @@
    "outputs": [],
    "source": [
     "number_steps = steps_to_gyroperiod * int(2 * np.pi)\n",
-    "trajectory = ParticleTracker(plasma, \"p\", 1, 1, timestep, number_steps)"
+    "trajectory = ParticleTracker(\n",
+    "    plasma, \"p\", 1, 1, timestep, number_steps, integrator=\"explicit_boris_relativistic\"\n",
+    ")"
    ]
   },
   {
@@ -210,6 +212,153 @@
     "trajectory.plot_time_trajectories()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can also take a look at the trajectory in the z direction:\n",
+    "    "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "trajectory.plot_time_trajectories(\"z\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "or plot the shape of the trajectory in 3D:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "trajectory.plot_trajectories()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As a test, we calculate the mean velocity in the z direction from the\n",
+    "velocity:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "vmean = trajectory.velocity_history[:, :, 2].mean()\n",
+    "print(\n",
+    "    f\"The calculated drift velocity is {vmean:.4f} to compare with the \"\n",
+    "    f\"expected E0/B0 = {-E0/B0:.4f}\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Relativistic variant"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "number_steps = steps_to_gyroperiod * int(2 * np.pi)\n",
+    "trajectory = ParticleTracker(\n",
+    "    plasma,\n",
+    "    \"p\",\n",
+    "    1,\n",
+    "    1,\n",
+    "    timestep / 100,\n",
+    "    number_steps * 50,\n",
+    "    integrator=\"explicit_boris_relativistic\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We still have to initialize the particle's velocity. We'll limit ourselves to\n",
+    "one in the x direction, parallel to the magnetic field B -\n",
+    "that way, it won't turn in the z direction.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "from astropy import constants\n",
+    "\n",
+    "trajectory.v[0][0] = 0.1 * u.m / u.s\n",
+    "trajectory.v[0][2] = 0.3 * constants.c"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Run the pusher and plot the trajectory versus time.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    },
+    "tags": [
+     "nbsphinx-thumbnail"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "trajectory.run()\n",
+    "trajectory.plot_time_trajectories()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

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
+
+
+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(