Adding a new physics package#

This guide will cover the files which need to be modified if you’re adding a new physics package. Small snippets of code will be included with the instructions to tell you where the new code needs to go. Once completed, you should be in a position to start writing functions which actually perform your new package.

physics_package.h#

Firstly, we will need a new alias for the physics package. The physics_package class will be used whichever package is active, but it will only call functions based on its type variable. We can add a new type by extending the enum class phys:

// Possible physics packages
enum class phys {
    bremsstrahlung,
    bethe_heitler_ee,
    bethe_heitler_muon,
    nuclear_trident,

Before we start the PIC loop, the physics package manager will call setup for each physics package in the list. This function switches to the setup function which matches the type variable. We need a corresponding setup function for our new package:

    void setup_bremsstrahlung(SimulationData<ndims> &data);
    void setup_bh_ee(SimulationData<ndims> &data);
    void setup_bh_mu(SimulationData<ndims> &data);

During the PIC loop, the physics package manager goes through each physics package, and calls the run command. Again, this switches to run the package which matches the type. We need to declare a run function for the new package in this file:

    // Run packages
    void run_bremsstrahlung(physics_package_manager<ndims> *manager,
                            SimulationData<ndims> &data);
    void run_bh_ee(physics_package_manager<ndims> *manager,
                   SimulationData<ndims> &data);
    void run_bh_mu(physics_package_manager<ndims> *manager,
                   SimulationData<ndims> &data);

Standard, non-collisional physics packages call two functions when they are run. A sampling function calculates the process rate, and determines whether it triggers for the current particle in the current time-step. A performing function actually modifies triggers macro-particles. We will need to declare functions for these, if appropriate

    void sample_bremsstrahlung(Species<ndims> *part,
                               bool &event_triggered,
                               physics_package_manager<ndims> *manager,
                               SimulationData<ndims> &data);
    void sample_bh_ee(Species<ndims> *part,
                      bool &event_triggered,
                      physics_package_manager<ndims> *manager,
                      SimulationData<ndims> &data);
    void sample_bh_mu(Species<ndims> *part,
                      bool &event_triggered,
                      physics_package_manager<ndims> *manager,
                      SimulationData<ndims> &data);

    void perform_bremsstrahlung(Species<ndims> *part,
                                physics_package_manager<ndims> *manager,
                                SimulationData<ndims> &);
    void perform_bh_ee(Species<ndims> *part,
                       physics_package_manager<ndims> *,
                       SimulationData<ndims> &);
    void perform_bh_mu(Species<ndims> *part,
                       physics_package_manager<ndims> *,
                       SimulationData<ndims> &);

If the process is collisional, then the sample/perform framework is inappropriate. Instead, sampling and performing must be done in the same function. The run function will loop over all cells, and a perform function must be declared for performing the process in the current cell

    void perform_coll_ir(int ix,
                         int iy,
                         int iz,
                         SimulationData<ndims> &,
                         physics_package_manager<ndims> *manager);

setup.cpp#

Now that we have a new physics package alias in the phys enum class, we need to tell the code when to use it. We need to add a new case to the setup_physics_packages function. This creates a new physics package for each physics_package block in the input deck. Here, the “val” string corresponds to the value given to the process key in the physics_package block. We need to link that input deck parameter to the type variable in the physics_package class

            if (val == "bremsstrahlung") {
                interaction->set_type(phys::bremsstrahlung);
            } else if (val == "bethe_heitler_ee") {
                interaction->set_type(phys::bethe_heitler_ee);
            } else if (val == "bethe_heitler_mu") {
                interaction->set_type(phys::bethe_heitler_muon);

physics_package.cpp#

Our physics_package class now has an alias for our new package. When the physics package manager calls a setup or run command for this class, we need to setup or run the corresponding package function. Hence, we need to add new cases to physics_package<ndims>::setup

    switch (type) {
        case phys::bremsstrahlung:
            setup_bremsstrahlung(data);
            break;
        case phys::bethe_heitler_ee:
            setup_bh_ee(data);
            break;
        case phys::bethe_heitler_muon:
            setup_bh_mu(data);
            break;

and physics_package<ndims>::run

    switch (type) {
        case phys::bremsstrahlung:
            run_bremsstrahlung(manager, data);
            break;
        case phys::bethe_heitler_ee:
            run_bh_ee(manager, data);
            break;
        case phys::bethe_heitler_muon:
            run_bh_mu(manager, data);
            break;

If your new package can create secondary particles, then there are additional functions which need to be modified. For each secondary species present, you will need to provide three pieces of information. Firstly, you must give the corresponding keyname for the species, as it appears in the input deck physics_package block. This is done in physics_package<ndims>::get_secondary_keynames

    switch (type) {
        case phys::bremsstrahlung:
            secondary_keynames.push_back("photon_species");
            break;
        case phys::bethe_heitler_ee:
            secondary_keynames.push_back("electron_species");
            secondary_keynames.push_back("positron_species");
            break;
        case phys::bethe_heitler_muon:
            secondary_keynames.push_back("muon_species");
            secondary_keynames.push_back("antimuon_species");
            break;

If the user doesn’t specify the secondary species in the input deck, the code will create a new species for this package. To do this, we need to provide the code with a default species name, given in the function physics_package<ndims>::get_new_spec_names:

    switch (type) {
        case phys::bremsstrahlung:
            new_species_names.push_back("bremsstrahlung_photon");
            break;
        case phys::bethe_heitler_ee:
            new_species_names.push_back("bh_electron");
            new_species_names.push_back("bh_positron");
            break;
        case phys::bethe_heitler_muon:
            new_species_names.push_back("bh_muon");
            new_species_names.push_back("bh_antimuon");
            break;

Finally, if a new species must be created, we must provide information about the type of species. A new case must therefore be added to physics_package<ndims>::get_new_spec_types.

        case phys::bremsstrahlung:
            new_species_types.push_back(part_type::photon);
            break;
        case phys::bethe_heitler_ee:
            new_species_types.push_back(part_type::electron);
            new_species_types.push_back(part_type::positron);
            break;
        case phys::bethe_heitler_muon:
            new_species_types.push_back(part_type::muon);
            new_species_types.push_back(part_type::antimuon);
            break;

Species.h#

Add any new phys_var variables needed. You should create a new optical depth variable to track when the new process triggers. The format currently in use is to add the od_ prefix for optical depth variables. If there are other permanent varibles needed, this is where to add them.

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

physics_package_manager.cpp#

Any new permanent particle variables we have added for the new package will need initialising. This is currently done in the function physics_package_manager<ndims>::initialise_shared_perm_vars. Make sure to add a new case for each new variable you have created.

    if (species_list[i]->od_bremsstrahlung.in_use()) {
        species_list[i]->od_bremsstrahlung.resize_array();
        for (auto part : *species_list[i]) {
            part->od_bremsstrahlung.set(-log(rng.random()));
        }
    }

For optical depth variables, we sample these from the emission probability distribution: \(\tau_d = -\ln\left(x_r\right)\), where \(x_r\) is a uniformly distributed random number. Note that this function will only initialise physics package variables which are “in use”. This is switched on by the packages themselves.

species.cpp#

Our new permanent variables need to be moved along with the rest of the particle data, whenever particle data is moved. Currently, this is done when macro-particles pass MPI boundaries, or when package-owned secondary particles are added to the main species list. Permanent physics package particle variables are packed for transfer in the function Species<ndims>::get_all_phys_variables, and we will need a new case for our new variables

    if (od_bremsstrahlung.in_use()) {
        phys_vars.push_back(od_bremsstrahlung.array_ptr());
    }
    if (od_bh_ee.in_use()) {
        phys_vars.push_back(od_bh_ee.array_ptr());
    }
    if (od_bh_mu.in_use()) {
        phys_vars.push_back(od_bh_mu.array_ptr());
    }

We also need to make that when packages create new particles, we initialise all the required permanent variables - including variables which needed for other processes. To do this, we pass each package-owned secondary species a list of all the active variables for the corresponding main-list species, using Species<ndims>::copy_phys_bools. Our newly added variables need to appear in this function too.

    if (other_spec->od_bremsstrahlung.in_use()) {
        od_bremsstrahlung.activate(*this);
    }
    if (other_spec->od_bh_ee.in_use()) {
        od_bh_ee.activate(*this);
    }
    if (other_spec->od_bh_mu.in_use()) {
        od_bh_mu.activate(*this);
    }

New package file#

For code neatness, it is good to create a new source-code file for each physics package. An easy way to implement this would be to find the test_process_*.cpp file which most closely resembles your new package, and to simply copy this file. Then, change the names of any setup, run, perform or sample functions to match those you defined in Species.h.

Once you reach this stage, you’ll have a new basic physics package which matches one of the test processes. From here, you can start modifying the test package functions to achieve the desired behaviour. This is how the real physics packages already implemented in the code were created. For details on more complicated implementations, try looking to the existing physics packages for guidance.

I’m aware elastic collisions are yet to be implemented. This may work differently to the current packages - look to the FORTRAN version of the code for help.

CMakeLists.txt#

All the source-code files appear in the SOURCES list of this file. For consistency, if you add a new package file, be sure to reference it here.

# Library sources exclude main.cpp, which contains main()
set(SOURCES
  Driver.cpp
  Fields.cpp
  Simulation.cpp
  Species.cpp

# Including headers in add_library calls, as this information is used by
# various IDEs.
set(HEADERS
  "${PROJECT_BINARY_DIR}/include/commit_info.h"
  "${PROJECT_SOURCE_DIR}/include/StaticArray1d.h"
  "${PROJECT_SOURCE_DIR}/include/StaticArray2d.h"
  "${PROJECT_SOURCE_DIR}/include/Array3d.h"

Adding a new physics package

Contents