NOTE: This section has not been updated in over a year and while much of the information is still accurate, some may be out of date. Please contact the developers with specific questions.

An electrostatic Particle-in-Cell (PIC) code for modeling Low Temperature Plasma devices. The code is designed from the ground up to be scalable and provides a base framework within which more advanced physics models can be incorporated.

Detailed Code Description

Table of Contents

Simulation Grid

Simulation Particles

Boundary Conditions

Communication

Source Code Files

C Structures & Memory Configuration

Input and Initialization

Time Loop

Simulation Grid

Table of Contents

The total simulation Grid is shown below. The grid is two-dimensional, uniform, structured and Cartesian with square cells. All grid based quantities are stored at cell nodes (also known as vertices). In all following figures, any area drawn on the lattice encompasses all of the nodes covered by its area, including those on its perimeter. The grid shown can be created by using the following settings under #Grid Properties within the input.dat file.

#Grid Properties
nx 42
ny 38
xmin 4.0
xmax 46.0
ymin 3.0
ymax 41.0

Note that this grid has 42 x 38 cells, which corresponds to 43 x 39 nodes.

If the input settings do not create square cells, the code will enforce this condition by altering ymax.

The Grid is broken up in to Regions (shown in orange), which are defined to handle particle decomposition in to different areas. Particles are communicated between Regions when they cross a Region boundary, however within a Region particles are parallelised via list splitting. That is, the list of particles is shared between all MPI tasks within a given Region.

Each Region is further broken up in to Domains (in blue). Each MPI task within the simulation owns a single Domain, which are defined to handle parallelisation for the field solve algorithm through Hypre.

Within each MPI task (aka within each Domain) there can be additional parallelisation via threading using OpenMP. A further level of parallelism is achieved via vectorisation.

The break down of the Grid in to Regions (orange) and subsequently Regions in to Domains can be created by the following settings under #Grid Properties within the input.dat file.

#Grid Properties
nx 42
ny 38
xmin 4.0
xmax 46.0
ymin 3.0
ymax 41.0
num_process_x 6 #Number of domains in the x-direction
num_process_y 4 #Number of domains in the y-direction
num_regions_x 2 #Number of regions in the x-direction
num_regions_y 2 #Number of regions in the y-direction

The only rules are that num_process_x % num_regions_x = 0 and num_process_y % num_regions_y = 0. That is the number of Domains (MPI tasks) in each direction must be divisible by the numer of Regions.

Boundary Conditions

Table of Contents

LTP-PIC is currently only capable of modeling a rectangular two-dimensional box. One can arbitrarily set four types of boundary conditions on each edge. These include:

• Periodic boundary conditions (option 0). This applies to both the field solver and simulation particles.
• Conducting boundary condition (option 1). This enforces a Dirichlet boundary condition for the field solver and an absorbing or exit boundary condition for the particles.
• Exit boundary condition (option 2). This enforces a Neumann boundary condition for the field solver and an absorbing or exit boundary condition for the particles.
• Mirror boundary condition (option 5). This enforces a Neumann boundary condition for the field solver and a reflective (mirror) boundary condition for the particles.

Periodic boundary conditions have additional restrictions which must be satisfied. These include that the opposite boundary must also be periodic, and that the number of cells in the periodic direction must be a power of 2. The code will warn you and exit if these conditions are not satisfied.

You can set the value of the electrostatic potential for the conducting boundary conditions via the input file (see the example below).

The boundary conditions for the grid used throughout this article can be set by the following settings within the input file:

#Boundary conditions, 0 - periodic, 1 - conducting, 2 - exit, 5 - mirror.
xbound_min 0
xbound_max 0
ybound_min 1
ybound_max 1
xbound_min_val 0.0 #Fixed potential for a conducting boundary condition.
xbound_max_val 0.0 #Fixed potential for a conducting boundary condition.
ybound_min_val 0.0 #Fixed potential for a conducting boundary condition.
ybound_max_val 1.0 #Fixed potential for a conducting boundary condition.

This sets up periodic boundary conditions in the x-direction (yes the power of 2 requirement has not been satisfied here - but I am too lazy to correct the figures), and conducting boundary conditions in the y-direction with a potential of 0 V at y-min and 1 V at y-max. Note that if a given boundary condition is not conducting, then the corresponding numerical value for the fixed potential is obviously ignored.

It is worth understanding how these boundary conditions effect the grid layout. Specifically, for this configuration, the Hypre field solver does not solve for the potential at the conducting boundaries (since this is known), nor does it solve for the potential at the right hand edge (since this this is identical to the value on the left hand edge). Therefore the grid which is constructed within Hypre covers the area shown in blue in the figure below, whereas particles can evolve over the entire grid (in red). Referring back to previous figures, you will see that the Domains are restricted to be only defined on this interior grid. Together the Domains cover the part of the grid which is fed to Hypre for the field solve. Note for exit and mirror boundary conditions, the edges are included in the field solve (and therefore are covered by the Domains).

Simulation Particles

Table of Contents

Work-in-progress

Communication

Table of Contents

Each Domain owns a collection of particles which can exist anywhere on its parent Region. However when a particle moves out of the Region boundary it must be communicated to a process in the new owner Region. Since the number and layout of Domains is identical for every Region, we have a logical one-to-one mapping between Domains on each Region. Particles are communicated between partner Domains on the different Regions as shown in the figure below.

To facilitate this communication we define a custom communicator called comm_part (see grid.c). The logical connectedness of this communicator is shown in the figure below. Note that this connectedness adheres to the grid boundary conditions such that our highlighted Domain is not connected to another process through the bottom, bottom left, or bottom right of its Region.

We also define a communicator to link all of the Domains on a given Region. This allows us to perform MPI Reduce operations within a Region. The communicator is called comm_region (see grid.c). The logical connectedness of this communicator is shown in the figure below.

We also have to consider the connectedness between each Domain and its 8 neighbours (regardless of their parent Region). This is necessary for the Hypre field solver, as well as sharing of edge data between Regions. A schematic of a Domain and its neighbours is shown in the following figure.

Source Code Files

Table of Contents

This is a list of the source files which make up LTP-PIC and a brief description of the functions they contain.

main.h

• Physical constants
• Code macro functions
• Structure definitions
• All functions headers

main.c

• Initialization of the simulation from input parameters.
• Main time loop.
• High level calls to Timing and diagnostics.

grid.c

• Initialization of the grid: Includes node layout, Regions, Domains, communicators, and memory allocation for grid based quantities.
• Checks for grid correctness.
• Charge summing.
• Edge communication routines for charge density, potential and electric field.
• Some grid based diagnostic routines.
• Clearing the grid.
• MPI_Reduce and MPI_Allreduce calls for grid quantities.
• Custom timed MPI_Barrier function.

field.c

• Initialization for Hypre. Setting up the grid, stencil and solver. Setting up charge density, potential and electric field vectors. Setting up Laplacian operator with boundary conditions and operators to compute electric field from the potential.
• Functions for mapping from Region grids to the Domains and visa versa.
• Poisson solver
• Computing electric field from potential.
• Calls to print Hypre operators and vectors.
• Finalizing Hypre

particles.c

• Memory allocation for each species and creation routine strucrures.
• Initialization of species, include allocating memory for particles.
• Initialization of creation routines, both of the plasma and injection type.
• Particle injection routines to be called at each time step.
• Particle memory handling routines.
• Particle push routines allowing sub-cycling and unmagnetized species.
• Particle boundary condition handling, including deletion, and communication between Regions.
• Sorting algorithms for particles.
• Interpolation routines for species density and electric field (zeroth and first order).
• Interpolation routines for diagnostics; particle flux and temperature (zeroth and first order).

numeric.c

• Initializing the Pseudo-Random-Number-Generator (PRNG) which utilised the Mersenne Twister algorithm.
• Random number generator for a Normal distribution.
• Other simple numerical functions used in grid setup.

dSFMT.c

• [External Package] functions for handling the Mersenne Twister PRNG.

timing.c

• Initialize Timing structure.
• Printing the time step during the loop.
• Accumulate final timing data and print to terminal and file.

input.c

• Defining analytic magnetic field.
• Reading the input and storing values in the respective strucures.
• Reading the species and creation input, and storing in the respective strucrures.

output.c

• Prepare directories and print initial simulation information.
• High level diagnostic function for printing particle phase data, density, current, temperature, charge density, potential, electric field
• Individual functions for printing each of the outputs.

checkpoint.c

• Routines for saving a checkpoint.
• Check if the simulation has been asked to start from a checkpoint.
• Allocate species memory and load from checkpoint.

landmark.c

• Routines specific to the Landmark 2a benchmark setup. Custom injection and handling of virtual cathode (not up to date).

C Structures & Memory Configuration

Table of Contents

Grid Structure

The grid structure, named Grid2D contains the following information:

• Logical grid layout, including the node limits owned by the entire Grid, by Regions, and Domains.
• Logical communication layout, including neighbouring Domains, neighbouring Regions and communicators.
• Grid boundary conditions at the edges of the entire grid, as well as Regions and Domains.
• Grid dimensions, and grid vectors giving the exact coordinates for the nodes (in SI units).
• Pointers to memory for stored grid quantities; for species density, charge density, potential, electric field, particle flux and temperature.
• Hypre containers for the grid, stencil, solver, operators and vectors.

Grid quantities which are species independent are accessed by a pointer to the allocated memory chunk. We set up our own logic for accessing grid components which sets the 0 index at the bottom left corner and cycles along the row (left to right) before moving up the column (bottom to top). The highest grid index is at the top right corner.

Grid quantities which are species dependent are accessed by a double pointer. This points to an array of arrays, the first depth of arrays has the size of the number of species. Each element of this array contains a grid memory chunk for storing the species property on the grid (density, flux, temperature etc.)

Each grid memory chunk is the size of the parent Region for the MPI process. No process owns the entire Grid (unless it is comprised of a single region), and only the Hypre vectors own memory chunks the size of the smaller Domain.

Species & Creation Structures

The AllSpecies structure contains information on the total number of species and a pointer to an array whose elements each contain a Species2D structure. Each Species2D structure contains the following information:

• Details on particle communication, including neighbouring regions.
• Species properties such as charge, mass and macro-particle weight.
• Modeling properties such as whether the particle is magnetized, or whether it undergoes sub-cycling.
• Local particle numbers and allocated memory for particles.
• Arrays for particle phase space data, electric field, cell identity and status.
• Basic grid information, including boundary conditions.
• Memory buffers for sending and receiving particles.

Data which is specific to each simulation particle are stored contiguously in arrays. That is, we rely on the Structure of Arrays (SoA) rather than Array of Structures (AoS) format. This should improve cache loading and vectorization performance.

The AllCreation structure contains the counts for the plasma creation and injection creation inputs and pointers to arrays of structures for CreatePlasma2D and CreateInjection2D.

The CreatePlasma2D structure contains information on:

• The species number.
• The number of particles to create.
• The area over which to create the particles.
• The initial particle temperature and directed kinetic energy.

The CreateInjection2D structure contains information on:

• The species number.
• The current density and number of particles to inject at each time step.
• The area over which to inject the particles.
• The injected particle temperatures and directed kinetic energy.

When initializing each Species2D, sufficient memory is created to store all of the particle properties for every registered CreatePlasma2D structure and an additional buffer is provided for CreateInjection2D structures. Each time injection occurs, a check is made on whether sufficient memory exists to store the new particles. If it does not exist, new memory is allocated with double the array length of the current allocation.

When particles leave a Region and are either deleted or communicated, a routine is run which replaces their particle information with that of the last particle in the list. This avoids memory holes at any time step and eliminates the need for logical conditions within particle push and interpolation routines.

The Species2D_buffer structure stores a subset of particle information for communicating between regions.

Other structures: Control, Timing, Output

The Control structure contains information on; the relative permittivity of the simulation space, the solver ID, restart logic and field sub-cycling.

The Timing structure contains information for recording timing profiles.

The Output structure contains flags for deciding which diagnostics to print during the diagnostic routine.

Input and Initialization

Table of Contents

The simulation begins by creating empty structures, counters and global timers. We then load and process the input file:

void ProcessingInput(struct Control *control, struct Grid2D *grid, struct AllSpecies *allspecies, struct AllCreation *allcreation, struct Timing *timer, struct Output *output);

This loads information on the time stepping control:

#Time Stepping Control
time_step_size 5.0e-12     #Time step size in s
number_steps 1000          #Number of time steps
print_interval 10          #Printing interval
sort_interval 100          #Interval over which to sort particles in memory
checkpoint_interval 100000 #Interval over which to store the simulation restart
field_subcycle 1           #Field sub-cycle interval

It also loads the grid properties, including processor layout and boundary conditions (see Simulation Grid for details):

#Grid Properties
nx 42           #Number of cells in the x-direction
ny 38           #Number of cells in the y-direction
xmin 4.0        #x-min in m
xmax 46.0       #x-max in m
ymin 3.0        #y-min in m
ymax 41.0       #y-max in m
num_process_x 6 #Number of domains in the x-direction
num_process_y 4 #Number of domains in the y-direction
num_regions_x 2 #Number of regions in the x-direction
num_regions_y 2 #Number of regions in the y-direction

#Boundary conditions, 0 - periodic, 1 - conducting, 2 - exit, 5 - mirror.
xbound_min 0       #x-min boundary condition
xbound_max 0       #x-max boundary condition
ybound_min 1       #y-min boundary condition
ybound_max 1       #y-max boundary condition
xbound_min_val 0.0 #Fixed potential for a conducting boundary condition.
xbound_max_val 0.0 #Fixed potential for a conducting boundary condition.
ybound_min_val 0.0 #Fixed potential for a conducting boundary condition.
ybound_max_val 1.0 #Fixed potential for a conducting boundary condition.

As well as details on models within the simulation:

#Models
relative_permittivity 1.0 #Relative permittivity over the entire simulation grid

It also stores information on which outputs the user might want to print:

#Output settings. 0 - do not print, 1 - print
print_particle_phase 0
print_densities 1
print_temperatures 1
print_currents 1
print_charge_densities 1
print_potential 1
print_fields 1

On this initial pass, the function records the total number of species and creation routines in order to prepare the appropriate amount of memory allocation.

Next we initialize the Mersenne Twister Pseudo-Random-Number-Generator (PRNG) seeding from /dev/urandom. A different seed is created for each MPI task and each OpenMP thread.

void InitializePRNG(dsfmt_t **dsfmt_handle);

Next we initialize all of the timers, setting them to zero.

void InitializeTimer(struct Timing *timer);

Next we check for a restart flag: ./pic -restart.

void CheckForRestart(int *argc, char **argv, struct Control *control, struct Timing *timer);

Next we initialize the grid, filling in all of the parameters for the Grid2D structure including the grid node placement, communication logic and memory for grid based quantities.

void InitializeGrid2D(struct Grid2D *grid);

Next we initialize the Hypre package. Setting up the grid, stencil, solver, operators and empty vectors.

int InitializeHypre(int argc, char **argv, struct Control *control, struct Grid2D *grid);

Next we allocate sufficient memory with AllSpecies to hold all of the Species2D structures, and sufficient memory in AllCreation for all of the CreatePlasma2D and CreateInjection2D structures.

void AllocateAllSpecies2D(struct AllSpecies *allspecies);

void AllocateAllCreation2D(struct AllCreation *creation);

Next we process and store the species input into each of the Species2D structures:

void ProcessingSpeciesInput(struct AllSpecies *allspecies);

An example of such an input and the input descriptions is given here:

#Species
species 1         #Species numeric identifier
clumping 1.28e6   #Clumping parameter, or macroparticle weight
charge -1.0       #Relative charge compared to the fundamental charge
mass 1.0          #Relative mass compared to that of an electron
magnetized 1      #Logic for whether the particle is magnetized or not
subcycle 1        #Sub-cycling rate for this species.
print_frac 0.0001 #Fraction of particles to print in the Phase Space diagnostic for this species


species 2
clumping 1.28e6
charge 1.0
mass 241097.065
magnetized 0
subcycle 1
print_frac 0.0001

Next we process and store the creation methods into the CreatePlasma2D and CreateInjection2D structures.

void ProcessingCreationInput(struct AllCreation *allcreation);

An example of such an input and the input descriptions is given here:

#Creation
plasma 1                #Plasma creation number
plasma_species 1        #Species number
plasma_particles 150000 #Number of macro-particles to create
plasma_xmin 0.0         #x-min for creation
plasma_xmax 0.0016      #x-max for creation
plasma_ymin 0.0         #y-min for creation
plasma_ymax 0.0024      #y-max for creation
plasma_temperature 10.0 #Initial species temperature in eV
plasma_kinetic_x 0.0    #Initial directed particle energy in the x-direction in eV
plasma_kinetic_y 0.0    #Initial directed particle energy in the y-direction in eV
plasma_kinetic_z 0.0    #Initial directed particle energy in the z-direction in eV

injection 1                      #Injection creation number
injection_species 1              #Species number
plasma_xmin 0.0                  #x-min for creation
plasma_xmax 0.0016               #x-max for creation
plasma_ymin 0.0                  #y-min for creation
plasma_ymax 0.0024               #y-max for creation
injection_temperature 10.0       #Initial species temperature in eV
injection_average_current -100.0 #Average injection current density (A/m^2)
plasma_kinetic_x 0.0             #Initial directed particle energy in the x-direction in eV
plasma_kinetic_y 0.0             #Initial directed particle energy in the y-direction in eV
plasma_kinetic_z 0.0             #Initial directed particle energy in the z-direction in eV

Next we initialize all of the species, allocating sufficient memory for particle storage. Then we proceed to particle creation including initial plasma species, and injection routines for creating future particles. If a restart has been requested, the particles are instead loaded from the Checkpoint directories.

void InitializeAllSpecies2D(struct Control *control, struct Grid2D *grid, struct AllSpecies *allspecies, struct AllCreation *creation, struct Timing *timer, dsfmt_t *dsfmt);

void InitializeAllCreation2D(struct Control *control, struct AllSpecies *allspecies, struct AllCreation *creation, struct Timing *timer, dsfmt_t *dsfmt);

Finally we create the Output and Checkpoint directories for all data to be saved during the simulation.

void PrepareOutputDirectories(struct Control *control, struct Grid2D *grid, struct Timing *timer);

void PrepareCheckpointDirectories(struct Control *control, struct Timing *timer);

We then proceed to the time loop!

Time Loop

Table of Contents

For describing the time loop we will mostly consider a single Domain which owns particles that can exist anywhere on its entire parent Region. Each Domain owns a copy of the grid quantities on the parent Region. Recall that this Domain is connected to its neighbouring Domains via mpi_comm_world, to the other Domains within its Region via comm_region and to its corresponding Domains within other regions via comm_part.

Begin by considering the collection of particles owned by our Domain, these particles could be of multiple species, although we visualise only one in our figure below.

Interpolating Particle Density to Grid

void InterpolateAllDensityToGrid2D(struct AllSpecies *allspecies, struct Grid2D *grid, struct Timing *timer);

In this first routine we interpolate the densities of the particles owned by our Domain on to the Region grid. The order of interpolation is set via compiler commands, with options:

• -DINTERP_ZERO - nearest grid point interpolation.
• -DINTERP_ONE - bilinear interpolation.

The function requires the AllSpecies, Grid and Timing structures. The following figure shows interpolation to the grid (particles from green to orange).

Summing Charge Density

void SumAllCharge2D(struct Grid2D *grid, struct AllSpecies *allspecies, struct Timing *timer);

Next we sum the total charge on the Region grid owned by our Domain. This loops through the grid for every species and adds their charge together.

The function requires the AllSpecies, Grid and Timing structures. The following figure shows the charge summing (particles from orange to red, and appearance of new particles from other species).

Collecting the Charge Density on the Region

void CommunicateGridChargeAllReduce(struct Grid2D *grid, struct Timing *timer);

Now we need to collect all of the charge on each Region by performing an MPI_Allreduce. This will sum all of the Region charge grids owned by each Domain and then send the total Region charge grid to each Domain (so each MPI task owns the complete copy).

This function requires the Grid and Timing structures and relies on the comm_region communicator (the connectedness of which is shown in the figure below).

Sharing Charge Edges

void CommunicateRegionGhostCellsCharge2D(struct Grid2D *grid, struct Timing *timer);

Now you may have noticed an issue. Since the Regions overlap at their edges some charge may have been interpolated to the same node owned by different Regions. We therefore need to share the Region edges. Studying the grid structure described earlier, we see that we only need to share and sum charge between the bottom to top and from the left to the right Region edges. To minimise cost, this is performed simultaneously by each MPI task, only sharing the part of the charge grid corresponding to its neighbour Domain across the Region boundary.

The function requires the Grid and Timing structures and relies on the mpi_comm_world communicator. The figure below illustrates this; charge is shared on the nodes surrounded by the dashed lines, and arrows show the direction that information is passed. Corners are automatically taken care of by ordering the sharing carefully.

Mapping the Region Charge to Hypre

int MapChargeToHypre2D(struct Control *control, struct Grid2D *grid, struct Timing *timer);

Now we map the charge on the Region charge grid to the Domain charge grid, or vector, which is owned by Hypre. This function requires the Control, Grid and Timing structure and is illustrated in the figure below.

Solving Poisson's Equation Using Hypre

int SolvePoisson(struct Control *control, struct Grid2D *grid, struct Timing *timer);

We let Hypre work its magic. The default solver is their PFMG Geometric Multigrid algorithm. Profiling consistently shows that it is the highest performing algorithm from the Hypre suite for this kind of global solve. This function requires the Control, Grid and Timing structure and requires the mpi_comm_world communicator (the connectedness of which is illustrated below).

Calculating the Electric Field Using Hypre

int ElectricFieldFromPotentialHypre(struct Control *control, struct Grid2D *grid, struct Timing *timer);

We similarly compute the electric field components from the potential via matrix multiplication in Hypre. This function requires the Control, Grid and Timing structure and uses the mpi_comm_world communicator (the connectedness of which is illustrated below).

Mapping the Domain Electric Field on Hypre to the Grid Regions

int GetFieldSolutionRegion(struct Grid2D *grid, struct Timing *timer); 

Now we perform the reverse operation as we did for the charge. We map from the Hypre vector (which owns the grid of a single Domain) to the Region grid owned by that task. Note that this will only map one Domain to the Region, we still need to collect all of these mappings in order to get the complete electric field on the Region grid.

This function requires the Grid and Timing structure and is illustrated in the figure below.

Sharing Electric Field Edges

void CommunicateRegionGhostCellsElectricField2D(struct Grid2D *grid, struct Timing *timer);

In order to interpolate back to the particles, we require the entire electric field on the Region grid. This means that the top and right edges of our Region will be lacking field information, since that was stored on the Domains from the corresponding Region on the other side of these edges. We share the electric field information from the top to the bottom and the right to the left at the Region edges. To minimise cost, this is performed simultaneously by each MPI task, only sharing the part of the electric field grid corresponding to its neighbour Domain across the Region boundary.

The function requires the Grid and Timing structures and relies on the mpi_comm_world communicator. The figure below illustrates this; electric field is shared on the nodes surrounded by the dashed lines, and arrows show the direction that information is passed. Corners are automatically taken care of by ordering the sharing carefully.

Collecting the Electric Field on the Region

void CommunicateGridElectricFieldAllReduce(struct Grid2D *grid, struct Timing *timer);

We now need to combine all of the electric field blocks which we have mapped to our Domain areas. This is again performed by an MPI_Allreduce, which will sum all of the electric field grids owned by each Domain and then send the combined Region electric field grid to each Domain (each MPI task).

The function requires the Grid and Timing structures and relies on the comm_region communicator (the connectedness of which is shown in the figure below).

Particle Injection

void AllInjection2D(struct AllSpecies *allspecies, struct AllCreation *creation, struct Timing *timer, dsfmt_t *dsfmt);

Additional particles are created if requested from an injection creation input. This function requires the AllSpecies, AllCreation and Timing structures, as well as a pointer to the PRNG. The figure shows the addition of new particles owned by the Domain in to the Region area.

Interpolate Electric Field to Particles

void InterpolateFieldToAllParticles2D(struct AllSpecies *allspecies, struct Grid2D *grid, struct Timing *timer); 

Now the electric field is interpolated back to the particles using the same order of interpolation set from the compiler flags. The function requires the AllSpecies, Grid and Timing structures. The following figure shows interpolation of field values to the particles (particles from orange to blue).

Particle Push

void PushAllParticles2D(struct AllSpecies *allspecies, struct Grid2D *grid, struct Timing *timer);

Particles are pushed by the standard Boris algorithm. Although the field solver is electrostatic, we can prescribe an analytic magnetic field. Heavier particles can also be pushed at intermittent time steps in order to reduce total computational cost. Furthermore we can use an unmagnetised field pusher if required.

The function requires the AllSpecies, Grid and Timing structures. The figure below shows the particles being moved slightly between time steps.

Particle Boundary Conditions

void HandleAllParticleBoundaries(struct AllSpecies *allspecies, struct Timing *timer);

Finally we handle particle boundary conditions for those particles which transition out of their parent Region during the time step. What happens to these particles depends on the boundary condition at the edge of the Region. For Periodic and Region to Region edges, the particles will be sent from their parent Domain to the corresponding Domain on the new Region (see the figure below). For conducting or exit Region edges, the particle will be deleted and removed from memory. For mirror Region edges the particle will be reflected.

The function requires the AllSpecies and Timing structures.

A typical time loop will then return to interpolation of the charge density back to the grid. However at certain time loops other routines are called. See below.

Particle Sorting

Diagnostics

Restarts

References