Physics package overview

Physics package overview#

The EPOCH physics packages are controlled using various different classes. We can see the relationships between classes at a glace, using the following diagram:

The main PIC loop is run by the Simulation<ndims>::solve function. The Simulation<ndims> class holds an instance of the physics_package_manager class, and this is called during the PIC loop. This is how the physics packages are accessed as the code runs.

physics_package_manager contains a physics_list variable, which is an array which holds physics_package classes. Every physics package set up by the input deck will feature in this list, and the code will run through each one individually at runtime.

Each physics_package is given an incident species ID, and the IDs of any background or secondary species. The physics_package class contains functions which can run any physics package, and a type variable is used to switch to functions which describe the current process. Originally, each package was a derived class of physics_package, but this approach was abandonned due to GPU concerns. For code clarity, functions related to each package are stored in different src files.

These physics packages require extra particle variables, which may be broken into 5 categories:

Number

Variable

Type

Stored in class

1

Temporary variable

One value per particle

Species

2

Permanent variable

One value per particle

Species

3

Species array

One value per cell

Species

4

Combined array

One value per cell

physics_package

5

Global array

One value per cell

physics_package_manager

  1. Temporary particle variables: sometimes multiple packages may need the same variable - e.g. total particle energy from momentum. We can calculate this once each step, and store this in the species class.

  2. Permanent particle variables: packages may need particle variables which are saved between steps, and are transferred with the particle when it moves between ranks. For example, most processes need an optical depth variable to track trigger points. This would also be needed for any “total work done by process” variable. A special phys_var class is used for these variables, which allows us to only allocate, move or calculate these variables if they are in use.

  3. Species background arrays: multiple packages may need a species number density. Save the full array if required to the species class.

  4. Combined-species background arrays: a physics package can have multiple background species, provided they all describe the same type of particle. These arrays sum the species background arrays, and save the result in physics_package.

  5. Global background arrays: some arrays require multiple different species to calculate, like the Coulomb logarithm. Arrays which can use the full species list are calculated and stored in the physics_package_manager.

The physics_package classes also have their own species list for storing newly created secondary particles. This prevents other processes from acting on new particles within the same step, so we can see their properties at creation in the output dumps. At the end of the step, these secondary particles are added back to the main list.

Physics package features#

The general capabilities of an EPOCH physics package are provided here:

All packages have a single incident species, which can interact with background fields, background species properties (like the number density), or in collisional processes where background macro-particles are affected.

When a package triggers, the incident particle may have its properties modified, or it may be destroyed. Secondary particles may be added to the package-owned species list, and background fields/macro-particles may be updated.

The trigger rate can be artificially increased or decreased in the physics package block. The process trigger rate can also be modified in a physical way, where the weights of secondary particles are increased or decreased to conserve the number of real secondary particles. For computationally intensive processes, we can improve code performance by subcycling, where we only run the package once every \(n\) steps, for a time-step \(n\) times longer. Additionally, calculation of the background arrays for packages can also be subcycled, if they are only expected to be slowly varying.

If the user doesn’t provide a secodary species for a package which needs one, the code will automatically generate a new species. Each package defines a default name and particle type for missing secondary species in the source-code.

Physics package manager class#

The physics package manager class, defined in physics_package_manager.h, is responsible for initialising and running each of the physics packages defined by the user in the input deck. A brief overview of the important features will be given here:

Temporary variables and arrays#

We create vectors in physics_package_manager<ndims>::setup_physics_packages to store a bitmask for each particle species:

   saved_inc_vars_per_species.resize(n_spec, 0);
   saved_bkg_vars_per_species.resize(n_spec, 0);
   remaining_inc_vars_per_species.resize(n_spec, 0);

Here, each element contains information about which temporary particle variables and arrays are needed for each process. For example, say saved_inc_vars_per_species[0] = 5. The binary representation of 5 is 101, which can be thought of as a logical array with the first and third elements set to true, and the second set to false. In this example, temporary variables [0] and [2] are active for species with ID 0, while [1] is not. The list of temporary variables is given in save_var.h:

constexpr std::array<var_pair, 13> save_var_map = {
    {{"ex", 1 << 0},
        {"ey", 1 << 1},
        {"ez", 1 << 2},

so saved_inc_vars_per_species[0] = 5 means ex and ez interpolated to the particle position are required for species 0, but not ey.

Physics packages cannot move the particle position, but the momentum can be changed. Variables derived from momentum (like energy part_E) need to be recalculated whenever the momentum is updated. However, if packages [0] and [1] need part_E, but no package after [1] uses it, then there is no point recalculating the temporary variable part_E if package [1] triggers. We recalculate the required-variable bitmask after each process triggers, and store this in remaining_inc_vars_per_species. This prevents us from updating temporary variables which are no longer required.

Global background arrays#

There are some global background arrays which could be accessed by multiple different packages, so we store these in the manager to allow each package access. Alternatively, the array may require species which are unknown to the package, as this only contains information for the incident species, and background species which must be equivalent to each other. Only the manager has access to all the particle species.

Currently, there are three global background arrays which can be calculated by the manager:

    // Arrays which hold global backgrounds
    Array3d<ndims> full_debye_length;
    void calc_debye_length(const std::vector<int> var_masks,
                           std::vector<Species<ndims> *> &spec_list,
                           const SimulationData<ndims> &data);
    Array3d<ndims> log_debye_length;
    void calc_log_debye_length(const std::vector<int> var_masks,
                               std::vector<Species<ndims> *> &spec_list,
                               const SimulationData<ndims> &data);
    Array3d<ndims> full_el_dens;
    void calc_el_dens(const std::vector<int> var_masks,
                      std::vector<Species<ndims> *> &spec_list,
                      const SimulationData<ndims> &data);

We require a calculator and an array variable for each one.

Permanent variables#

Individual physics packages activate the permanent particle variables they need, but the manager is the class which initialises active permanent variables. The `` function performs two tasks. Firstly, the active permanent variables in the main species list are also switched on for the corresponding species in package-owned lists:

    // Copy required permanent variables to package-owned species lists
    for (auto &interaction : physics_list) {
        interaction->secondaries_setup(data);
    }

This is required for particle transfer from package-owned species lists to the main species list. Secondly, the initial values for the permanent variables are set. These values can be anything, although most permanent variables currently in the code are optical depths, which are initialised from random sampling of the trigger optical depth.

    // Set-up and initialise permanent particle variables shared between physics
    // packages
    initialise_shared_perm_vars(data.species_list);

Particle binning#

In collisional processes, we allow modification of both incident and background macro-particles. To enforce local-collisions only, we bin particles by creating an array of vectors in the Species class Array3d<ndims, std::vector<int>> part_n_in_cells. Each element in the array gives a vector of particles in the current cell. In the Species class, the properties of different particles can be extracted by changing the integer n. Hence, part_n_in_cells gives a list of particle n values for the current cell.

Additional Species functions are provided for navigating through binned particles. The function bool cycle_through_cell(ix,iy,iz) sets the n value to the top of the list for cell (ix,iy,iz), or returns false if the cell is empty. The function next_part_in_cell() sets n to the next particle in the cell, or returns false if there are no more particles in the list. To cycle through binned particles, we must call cycle_through_cell, before calling next_part_in_cell.

Running packages#

When Simulation<ndims>::solve is running the PIC loop, we access the physics packages by running: physics_package_manager<ndims>::run_physics_packages. This function does quite a few things, so we will write a quick summary of the order of operations, along with the motivations.

  • Firstly, loop through all physics packages and determine which are active this step. The input deck allows users to run physics packages once every \(n\) steps.

  • Determine which temporary variables need to be pre-calculated this step, then calculate them.

  • For species involved in collisional processes, bin these to cells.

  • Calculate the global background arrays, if any package needs them.

  • Create combined-background arrays for processes with multiple equivalent background species.

  • Cycle through every package after the current one in the list, and figure out which temporary variables they need. This will determine which ones get re-calculated if the incident kinematic variables change.

  • Interpolate the global background arrays to the incident particle positions, if needed.

  • Run the package

Once all packages have run, this function copies particles from the package-owned secondary species to the main list.

Physics package class#

The physics_package class contains many functions, which can be broadly broken down into a few categories:

  • Getters and setters for input deck keys, like sampling parameters, and species ID for incident, background and secondary species

  • Functions to setup package-owned the species list, and to initialise and transfer new particles

  • Variables and functions related to combining background arrays

  • Package-specific setup/run functions

  • Package-specific variables/functions for rate calculation

  • Particle transformations: momentum rotation, Lorentz boost and inverse boost

Some of these features are quite straightforward, especially with the context of the physics packagage manager. This section will focus on a few specific points

Package-specific variables#

There are a few parameters which need to be specified in package-specific setup functions. For example, some of the lines in setup_bremsstrahlung read:

    inc_vars_to_save = 0;
    inc_vars_to_save |= var_code("g_facs");
    inc_vars_to_save |= var_code("part_v");
    inc_vars_to_save |= var_code("part_E");
    inc_vars_to_save |= var_code("part_ni");
    bkg_vars_to_save = 0;
    bkg_vars_to_save |= bkg_code("num_dens");

    auto inc_spec = data.species_list[incident_id];
    inc_spec->od_bremsstrahlung.activate(*inc_spec);

We must construct the bitmasks for inc_vars_to_save, and bkg_vars_to_save. We can do this using the bitwise-or command |=, along with the incident variable and background array variable codes found in save_var.h. Here, g_facs refers to the macro-particle interpolation factors on the grid, which are used to turn the background number density array num_dens into a number density at particle variable part_ni.

Additionally, we use this function to activate any permanent particle variables which are needed for the species. Activating a variable does not initialise it, but it marks it as a variable which needs initialising by the physics package manager.

Combined arrays#

In the FORTRAN code, if the user ran a bremsstrahlung simulation with an electron species passing through two different neutral gold species, then bremsstrahlung would be run twice for each background. We can speed up the code by combining equivalent species together.

In this case, the bkg_vars_to_save bitmask for "num_dens" would be set for both background species, calculated, and stored in the Species variable. Species-owned background arrays are stored in a Species map: std::unordered_map<int, Array3d<ndims>> background_arrays;, and can be accessed using the same integer key used in bkg_vars_to_save:

    if (var_mask & bkg_code("num_dens")) {
        Array3d<ndims> *arr = &species->background_arrays[bkg_code("num_dens")];
        arr->resize(-ng, i1, -ng, j1, -ng, k1);
        species->calc_num_density(*arr);
    }

This particular form was chosen as it means we only need to create map elements for arrays which are currently in use. However, I think map lookup may be slow, and we might want to move to creating something similar to phys_var for the permanent particle variables.

If multiple backgrounds are used, the package will combine all these background arrays together, and save them in a package-owned array:

    // Arrays which hold summed background properties
    Array3d<ndims> full_bkg_num_density;
    Array3d<ndims> full_bkg_temperature;

For single-background packages, the array in the Species map is used. For multi-background packages, the physics_package array is used. To transition between the two, a pointer is used to point to whichever array is appropriate for the current package.

    // Pointer to background arrays, either for full or single-species arrays
    Array3d<ndims> *bkg_num_density;
    Array3d<ndims> *bkg_temperature;

physics_package block#

The physics_package block in the input deck may be written in the format:

begin:physics_package
  process = bremsstrahlung
  incident = Positron_bunch
  background = Gold
  photon_species = Photon_from_pos
  min_secondary_energy = 0
end:physics_package

The block takes the following keys:

  • process: The name of the physics package process to use. This is used to set the package type.

  • incident: The name of the incident particle species.

  • background: The name of the background species. Multiple background species can be provided, as long as they’re all equivalent species (same mass, charge, atomic number).

  • Secondary species: The user can provide the names of secondary species. Different physics packages set different secondary particle types, so different keys are used here. Acceptable secondary species keys include: secondary1, secondary2, secondary3, photon_species, electron_species, positron_species, muon_species, antimuon_species, ejected_electron, ionise_to_species, recombine_to_species

  • incident_recoil: When a process triggers, modify the incident particle species if the package is capable of doing so. Default is T, which is the correct physical behaviour. This can be switched off for debugging/testing purposes.

  • produce_secondaries: If a process creates secondary particles, add the macro-particles to the simulation. Default is T, which always creates secondary macro-particles if the pacakge is capable of doing so. This can be switched off for debugging/testing purposes.

  • run_on_nth: Only run the process once every \(n\) steps, for a time-step of \(n\Delta t\). This key accepts an integer parameter. The default is 1, to run every step.

  • artificial_rate_boost: Artificially modify the process trigger rate by some factor. The default is 1, for no unphysical boost. Numbers greater than 1 will increase the process rate, and numbers less than 1 decrease the rate.

  • boost_secondaries: Modify the process trigger rate by some factor, but change the weights of secondary macro-particles to conserve the real particle number. The incident particle trigger response will be modified to ensure no change in the incident physical behaviour. The default is 1, for no weight change. Numbers greater than 1 yield more common secondary emission, with lower weights. Numbers less than 1 yield fewer secondary emissions, with greater weights.

  • recalc_time_num_density: Specifies how long we wait before recalculating the background number density. Default is 0, to recalculate every step.

  • recalc_time_temperature: Specifies how long we wait before recalculating the background temperature. Default is 0, to recalculate every step.

  • recalc_time_debye_length: Specifies how long we wait before recalculating the background Debye length. Default is 0, to recalculate every step.

  • recalc_time_log_debye: Specifies how long we wait before recalculating the natural logarithm of the Debye length. Default is 0, to recalculate every step.

  • total_pz_transfer: Parameter for test_process_EB_background. Maximum allowed increase of pz to a single particle.

  • min_secondary_energy: Secondary macro-particles are only added to the simulation if the sampled secondary energy is greater than this number.

  • laser_omega: Parameter for field ionisation, as these rates depend on the laser angular frequency.

  • use_brem_plasma_screening: Key to switch on a plasma-screening correction to the bremsstrahlung process.

  • sample_angular_dist: If active, some packages will use a more complex model for the angular distribution of emitted secondaries.

  • use_collision_ionise: If using the collisional ionisation and recombination package, this switches on/off the collisional ionisation part. Default T.

  • use_three_body_recombine: If using the collisional ionisation and recombination package, this switches on/off the three-body recombination part. Default T.

  • use_dielectronic_recombine: If using the collisional ionisation and recombination package, this switches on/off the dielectronic recombination part. Default T.

  • use_radiative_recombine: If using the collisional ionisation and recombination package, this switches on/off the radiative recombination part. Default T.

  • use_multiphoton: If using field ionisation, this switches on/off the multiphoton ionisation part. Default T.

  • use_tunnelling: If using field ionisation, this switches on/off the tunnelling (ADK) ionisation part. Default T.

  • use_barrier_suppression: If using field ionisation, this switches on/off the barrier-suppression ionisation part, using the classical Posthumus method. Default T.

  • use_classical_rate: If using field ionisation, this switches on/off the barrier-suppression ionisation part, using the simple classical method. Default is F. If you want to use this, you also have to set use_barrier_suppression = F.

Physics table class#

For computational efficiency, rates and other parameters related to sampling physics packages are pre-calculated in tables, and interpolated values are returned. The physics_table class in physics_table.h was designed to provide a general interface for the tables used in all physics packages.

Initially, the format was quite straightforward. The table contained general 1D and 2D array parameters:

    unsigned int nx;
    unsigned int nz;
    std::vector<real> x_vals;
    std::vector<real> y_vals;
    std::vector<std::vector<real>> z1_vals;
    std::vector<std::vector<real>> z2_vals;
    table_type type;

and different packages would have different functions to initialise relevant arrays using package-specific tables:

    void read_brems_table(const std::string &phys_path, int atomic_no);
    void setup_bh_ee_tables(real bh_screen_fac, real bh_fz);
    void setup_nuc_trident_tables();
    void setup_ci_tables_mbell(std::vector<real> bind_ene,
                               std::vector<real> el_n,
                               std::vector<real> el_l,
                               int atomic_no);

Once the data was stored in x_vals, y_vals and the 2D z1_vals, z2_vals, the code could then use generic interpolation functions. For example, in the function:

real get_y_from_x_y(const real &x_in);

we assume \(x\) and \(y\) are equally-long 1D arrays, and we perform a simple interpolation to get \(y(x_\text{in})\). 2D interpolators were also given, which act by performing a series of multiple 1D interpolations.

As more and more packages were added, this simple setup was no longer appropriate. Some tables have arbitrarily-spaced sample points, while some have logarithmically-spaced points. The arbitrary tables need bi-section interpolation, while direct look-up is faster for the logarithmic ones. Some packages also need many tables. As a result, this class has become quite messy. I think it would be better if each package had its own derived table class, but I talk about this more in “unfinished work”.

Particle variable class#

In the Species class, we define a new phys_var variable each time we want to add a new permanent particle variable. Here, the phys_var class solves a memory issue. Particles in the Species class are represented by arrays, where Species has an internal variable n, and particle access functions return the nth element of the array. But when we start adding dozens of new particle variables for the physics packages, allocating these arrays will start to hurt the memory.

Hence, permanent variables are given the phys_var type. These have a boolean key to determine whether the variable is needed or not. If it is needed, then memory will be used to create an array for the current species, holding the particle data for this variable. It is defined as a friend class in Species, such that it can access the current n:

    // Permanent MPI-aware particle variables
    friend class phys_var<ndims>;
    phys_var<ndims> od_bremsstrahlung;
    phys_var<ndims> od_bh_ee;
    phys_var<ndims> od_bh_mu;
Contents