Added relativistic boris push

[pheuer]

Added relativistic Boris push to simulation.particle_integrators. This is necessary for the synthetic proton radiography module (protons are typically > 1% c).

I think I wrote it correctly (based on Birdsall), but it would be great if someone could check my work.

Also, currently there are no tests. Ideally this should be tested as part of the particletracker tests (maybe showing that it produces the same result in the non-relativistic limit as the regular Boris push) but those tests are currently in flux.

[StanczakDominik]

A few suggestions, and we do have to figure out how to test this, because the implementation looks correct enough to fool me 😄 I did start reworking the particletracker tests in #675, and I'll prioritize bringing them in for now. To get this merged before I do, I see two options atm:

1. Hack together some tests like, "two trajectories basically equivalent at low initial velocities, to floating point accuracy", "two trajectories completely different according to some heuristic such as end position/velocity at high initial velocities"
2. Wait until we get some better ideas and move this into _wip_integrators in particletracker; also add a warning at https://github.com/pheuer/PlasmaPy/blob/relativistic_boris_push/plasmapy/simulation/particletracker.py#L113 if integrator not in self.integrators that this one is untested for now. Possibly add a filterwarnings in your code to ignore that.

Whatcha think?

[StanczakDominik]

Note to self: look into this!

[pheuer]

Yes please do - I'm very skeptical whenever my math disagrees with the textbook!

[pheuer]

@StanczakDominik I like option 2: saves us from writing tests now that will just be replaced with better ones soon. I've added a warning like you suggested - what do you think?

[StanczakDominik]

Sounds good! I'll dig into that paper as soon as I'm able, and we can probably merge it.

[pheuer]

@StanczakDominik Thanks - I think I was just looking at the definitions in section 2.2.1 of that paper (although maybe the Lorentz-invariant formulation they talk about later is worth looking into).

[StanczakDominik]

Okay, change of plans due to telecon where I was reaffirmed in thinking that having this go in untested is a meh idea; I'll add them.

• look through references, find possible test cases
• copy tests for other integrator
• symplectic integrator so should be reversible?

[StanczakDominik]

And if you wanna continue building on top of this, just do so and we'll do git magic to account for it :)

[StanczakDominik]

Here's a screenshot from the talk I watched, that I mentioned yesterday at the telecon:

[codecov]

Codecov Report

Merging #979 (a22e24e) into master (af7b9d3) will decrease coverage by 0.79%.
The diff coverage is 26.31%.

@@            Coverage Diff             @@
##           master     #979      +/-   ##
==========================================
- Coverage   96.86%   96.07%   -0.80%     
==========================================
  Files          70       70              
  Lines        6864     6881      +17     
==========================================
- Hits         6649     6611      -38     
- Misses        215      270      +55     

Impacted Files	Coverage Δ
plasmapy/simulation/particle_integrators.py	53.84% <14.28%> (-46.16%)	⬇️
plasmapy/simulation/particletracker.py	36.06% <60.00%> (-62.22%)	⬇️
plasmapy/plasma/sources/plasma3d.py	93.75% <0.00%> (-6.25%)	⬇️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update af7b9d3...a22e24e. Read the comment docs.

[StanczakDominik]

So this is pretty wack. I'm pretty sure I said somewhere that boris_push takes ndarrays, not Quantities.

Looks like I lied! The tests are failing now, and the way I've got them fixed locally is actually using _c = constants.c.

This means we've got a pretty simple speed optimization available before us; but it's out of scope of this PR, so let's get this one merged faster. I'll continue reworking the particle tracker and optimize that out.

[StanczakDominik]

I just want to restate how confused I am right now.

[StanczakDominik]

Okay, I figured it might be a good idea to still test whether the relativistic case outputs correct results.

Doesn't look like it, based on the notebook.

And since the particletracker code is a bit of a dumpster fire right now (with myself, of course, to git blame), I think it makes more sense to come back to this once I've got that cleaned up.

[StanczakDominik]

Well, admittedly I don't even know good test cases for the relativistic pusher...

[pheuer]

I agree, we can come back to this a bit later.

What test did you run on the relativistic pusher (the notebook you mentioned earlier)? One option would of course be making sure that it matches the non-relativistic Boris push for non-relativistic particles. Another might be pushing a relativistic particle a single time step in a known uniform field and somehow analytically calculating what the right answer should be another way to compare?

[StanczakDominik]

I did run the non-relativistic limit as a test, and the results agree; but I'd feel bad having a relativistic pusher go in without testing it in the relativistic limit...

And I have little subject matter expertise on the relativistic limit, so, yeah.

[pheuer]

@StanczakDominik So I'm not sure how to best integrate this into test_particletracker, but I've come up with a good test problem. This very nice paper gives the full solutions for ExB drift of a relativistic charged particle. Currently I've only implemented the case where E<cB, but the paper has E=cB and E>cB too. The paper also includes the analytical non-relativistic solution.

In the following plots I'm plotting the coordinate in the direction of the drift (parallel to E). If E=0.05 * cB, the ExB drift is effectively non-relativistic:

But if E = 0.5 * cB, there is a clear relativistic effect

I'm not entirely sure why vd=E/B matches the relativistic prediction here...I thought that was the classical ExB drift?

Would it be difficult to implement a test that runs particle_tracker with explicit_boris_realtivistic and checks the output trajectory against this analytical solution?

Here's the code:

import astropy.units as u
import astropy.constants as const
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

# This function is Eq. 72 for case 1
def t_opt(𝜏, t, vd, γd, ν):
    """
    fsolve optimization function
    Eq. 72 in the Friedman paper

    """
    return t - γd.value**2/ν.value*(ν.value*𝜏 -  vd.value**2/const.c.si.value**2*np.sin(ν.value*𝜏))

def ExB_trajectory(t, E, B, q=const.e.si, m=const.m_p.si,
                   nonrel=False):
    """
    Calculates the relativistically-correct ExB drift trajectory for a
    particle starting at x=v=0 at t=0 with E in the x direction and
    B in the y direction and ExB in the z direction. The motion is 2D in the
    xz plane.

    From https://journals.aps.org/pre/abstract/10.1103/PhysRevE.72.026603

    """

    # TODO: implement this for more regimes?
    if E >= const.c.si*B:
        raise ValueError("Currently this function only works for E<cB")

    # Just after Eq. 41
    vd = (E/B).to(u.m/u.s)

    # Eq. 34
    ν = q*np.sqrt((const.c.si*B)**2 - E**2)/(m*const.c.si)

    if nonrel:
        # Calculate the positions using the non-relativistic expression
        # Eq. 79
        a = m*vd/(q*B)
        x = (a * (ν*t - np.sin(ν.si.value*t.si.value))).to(u.m)
        z = (a * (np.cos(ν.si.value*t.si.value) - 1)).to(u.m)
    else:
        # Eq. 45
        γd = 1/np.sqrt(1 - vd**2/const.c.si**2)

        # Numerically invert Eq. 72 to calculate the proper time for each time
        # t in the lab frame
        𝜏 = np.zeros(t.size)
        for i, val in enumerate(t):
            𝜏[i] = fsolve(t_opt, [val.si.value], args=(val.si.value, vd, γd, ν))
        𝜏 *= u.s

        # Calculate the positions
        # Eq. 71
        a = γd*vd/ν
        x = (a*γd*(ν*𝜏 - np.sin(ν.si.value*𝜏.si.value))).to(u.m)
        z = (a*(np.cos(ν.si.value*𝜏.si.value) - 1)).to(u.m)

    return x, z



# Timescale is ~ns because fci is ~6 ns
t = np.linspace(0, 100, num=100) * u.ns

# Three regimes: E<cB, E=cB, E>cB
regime = 0.5

B = 10 * u.T
E = (regime * const.c.si * B).to(u.V/u.m)



x_rel, z_rel = ExB_trajectory(t, E, B)
x_nonrel, z_nonrel = ExB_trajectory(t, E, B, nonrel=True)

# Calculate the simple vd = E/B position predictions
x_nonrel_linear = ((E/B) * t).to(u.m)

fig, ax = plt.subplots()

ax.plot(t, x_nonrel_linear, label='vd=E/B')
ax.plot(t, x_nonrel, label='Non-relativistic')
ax.plot(t, x_rel, label='Relativistic')
ax.set_xlabel("t (ns)")
ax.set_ylabel("x (m)")
ax.set_title(f"E/cB = {regime}")
ax.legend()

[StanczakDominik]

Oh hell yeah, that's a great test! It shouldn't be too hard; stick an L2 norm on the difference between the two trajectories and it should (hopefully...) disappear (or at least, become very smol) for small timesteps. I can try handling that after Friday.