Unfinished work

Unfinished work#

Packages#

I was unfortunately unable to port over the following packages in time:

  • Collisions (within 1 species)

  • Collisions (between 2 species)

  • Background collisions (between 2 species, only affecting an incident species)

  • Non-linear Compton scatter (qed)

  • Breit-Wheeler pair production (qed)

  • QED trident pair production

The collisions packages are unlike any of the other packages currently present in the code. Since the time between collisions is much smaller than the EPOCH time-step, we don’t sample each individual collision. Instead, the code estimates how many collisions are performed between calls, and samples a scatter angle based on this. When you port over collisions from the FORTRAN code, you might want to look into the correction described in Higginson (2020), “A corrected method for Coulomb scattering in arbitrarily weighted particle-in-cell plasma simulations”. Journal of Computational Physics, 413, 109450.

Nonlinear Compton scatter follows the same framework as bremsstrahlung radiation, Breit-Wheeler follows Bethe-Heitler, and the two trident packages are similar. The only difference is that the rates of these processes depend on the background field, not background ion densities. I put these packages last, as I’m aware there are efforts to update these in York first.

Multiple triggers per timestep#

Most packages will only trigger once every few time-steps, so an ability to allow multiple triggers per time-step isn’t required. Collisions trigger so many times within a time-step that it isn’t practical to model each one individually

  • this is why we use a cumulative scatter angle.

However, there is one package which triggers too often to ignore multiple-triggers, but not often enough to switch to a cumulative trigger option: field ionisation. In the FORTRAN code, an ion can ionise multiple times within a single timestep - this is the only package which allows multiple sampling.

One way we could implement this would be as follows:

  • For each package, make a new frac_dt variable - this tracks the fraction of the time-step sampled. Initially, frac_dt = 1 for every particle, every step.

  • Perform all physics packages. If a particle doesn’t trigger, frac_dt goes to 0. If it does trigger, then the process optical_depth has gone negative.

  • If a package triggers, set frac_dt to ABS(optical_depth) / delta_optical_depth. For example, if an optical depth goes from 0.5 to -1.0 over a time-step, we say that the interaction happened after one third of a time-step, with two-thirds of a time-step left to sample

  • Move all package-owned secondary particles to the main particle list - these have a frac_dt which match the incident particle. Note that some incident particles will be destroyed at this stage - set secondary frac_dt to what the incident would have been.

  • Perform all physics packages again, but only for particles with frac_dt > 0

This method has the advantage of being applied to every physics-package, so we’ll have the correct behaviour not just for field ionisation, but for any process which triggers multiple times per timestep. This relaxes the constraint on upsampling rare processes, as we don’t have to worry as much about hitting the one event per timestep limit. On the downside, we could be trading this accuracy improvement for a performance hit. Maybe allow the user to switch off the multiple-trigger check for packages in the physics_package block.

If we did implement this, we would need to change the optical depth sampling behaviour. Currently, we add a newly sampled optical depth to the old negative optical depth value when a process triggers. In this way, the extra optical depth travelled in the previous step is used as distance towards the next trigger. With this new method, that extra optical-depth is captured by the frac_dt variable - so that we don’t double-count this, we would need to set the optical depth to 0, and then add the newly sampled optical depth.

Master blocks#

In the FORTRAN code, there was no physics_package block. Users would set physics packages by using things like a bremsstrahlung block, or a qed block. To make the input decks backwards-compatible, we would need to read these old master blocks, and use them to automatically generate physics_package blocks.

For example, the bremsstrahlung block would look for every electron species with identify:electron (or equivalent electron identity), and mark these as incident species. The code should then scan for every background species marked with an atomic number. Once we have a list of electron and background species, the code should create a unique physics_package instance for each incident/background combination, and add all of these to the physics list. There was no consistent design-style for these master blocks, so each would need to be implemented separately.

The purpose of these master blocks would be to save time setting up each physics_package individually. However, if the user did manually specify each physics package, then they would have a greater control over which ones are sub-cycled, or they could have different secondary species populated for different interactions. There is greater flexibility without the maser-blocks, so I didn’t consider them a priority for my development.

Output/restart#

The physics packages introduce new particle variables, some of which need to be saved between steps. Currently, there is no way to output these new variables, which means that restart dumps would not be possible. Parameters like “optical depth” are required for restarts.

Python interface#

Currently, physics packages can only be set up from the input deck - there is no support for the Python interface.

Tidy up tables#

I think the structure of the physics_table class could be improved. The idea behind this was to structure our data so we only ever needed to implement one table interpolator function, for each type of interpolation. However, different packages have different numbers of tables, and we end up with a mess of functions like get_z3_from_log_y_z4, which makes no sense on its own.

I actually think it would be better if there was a different table class for each package. These would be stored in the physics_package class, but only used if the package was the right type. Then we could use sensible names for each of the tables present, and these classes could call functions from a separate class which only contained interpolation functions like interpolate_bisection.

Contents