Write a namelist¶
Before you run Smilei, you need a namelist (an input file). The namelist is written in the python language. It is thus recommended to know the basics of python.
We suggest you copy one existing namelist from the folder benchmarks.
All namelists have the extension .py
.
General rules¶
Smilei requires a few blocks to be defined, such as:
Main( # ... timestep = 0.01, # defines the timestep value grid_length = [10., 20.], # defines the 2D box dimensions # ... )
Outside blocks, you can calculate anything you require. Inside a block, you must only define variables for Smilei.
The python syntax requires special indentation of each line. You begin with no indentation, but you have to add four spaces at the beginning of lines inside a group, and so on. For instance:
if a == 0: timestep = 0.1 if b == 1: timestep = 0.2 else: timestep = 0.3
You will need to use lists, which are series of things in python, defined between brackets
[]
and separated by commas. For example,mean_velocity = [0., 1.1, 3.]
.You are free to import any installed python package into the namelist. For instance, you may obtain \(\pi\) using
from math import pi
.All quantities are normalized to arbitrary values: see Units.
Python workflow¶
Python is started at the beginning of the simulation (one python interpreter for each MPI process). The following steps are executed:
A few variables from Smilei are passed to python so that they are available to the user:
The rank of the current MPI process as
smilei_mpi_rank
.The total number of MPI processes as
smilei_mpi_size
.The maximum random integer as
smilei_rand_max
.
The namelist(s) is executed.
Python runs
preprocess()
if the user has defined it. This is a good place to calculate things that are not needed for postprocessing with happi.The simulation is initialized (including field and particle arrays).
Python runs
cleanup()
if the user has defined it. This is a good place to delete unused heavy variables.Python checks whether the python interpreter is needed during the simulation (e.g. the user has defined a temporal profile which requires python to calculate it every timestep). Otherwise, python is stopped.
All these instructions are summarized in a file smilei.py
,
so that the user can directly run python i smilei.py
for postprocessing purposes.
Main variables¶
The block Main
is mandatory and has the following syntax:
Main(
geometry = "1Dcartesian",
interpolation_order = 2,
interpolator = "momentumconserving",
grid_length = [16. ],
cell_length = [0.01],
simulation_time = 15.,
timestep = 0.005,
number_of_patches = [64],
cluster_width = 5,
maxwell_solver = 'Yee',
EM_boundary_conditions = [
["silvermuller", "silvermuller"],
# ["silvermuller", "silvermuller"],
# ["silvermuller", "silvermuller"],
],
time_fields_frozen = 0.,
reference_angular_frequency_SI = 0.,
print_every = 100,
random_seed = 0,
)
 geometry¶
The geometry of the simulation:
"1Dcartesian"
"2Dcartesian"
"3Dcartesian"
"AMcylindrical"
: cylindrical geometry with Azimuthal modes decomposition.
In the following documentation, all references to dimensions or coordinates depend on the
geometry
. 1D, 2D and 3D stand for 1dimensional, 2dimensional and 3dimensional cartesian geometries, respectively. All coordinates are ordered as \((x)\), \((x,y)\) or \((x,y,z)\). In the"AMcylindrical"
case, all grid coordinates are 2dimensional \((x,r)\), while particle coordinates (in Species) are expressed in the 3dimensional Cartesian frame \((x,y,z)\).Warning
The
"AMcylindrical"
geometry has some restrictions. Boundary conditions must be set to"remove"
for particles,"silvermuller"
for longitudinal EM boundaries and"buneman"
for transverse EM boundaries. You can alternatively use"PML"
for all boundaries. Vectorization, collisions and order4 interpolation are not supported yet.
 interpolation_order¶
 Default
2
Interpolation order, defines particle shape function:
2
: 3 points stencil, supported in all configurations.4
: 5 points stencil, not supported in vectorized 2D geometry.
 interpolator¶
 Default
"momentumconserving"
"momentumconserving"
"wt"
The interpolation scheme to be used in the simulation.
"wt"
is for the timestep dependent field interpolation scheme described in this paper .
 grid_length¶
 number_of_cells¶
 A list of numbers: size of the simulation box for each dimension of the simulation.
Either
grid_length
, the simulation length in each direction in units of \(L_r\),or
number_of_cells
, the number of cells in each direction.
 cell_length¶
A list of floats: sizes of one cell in each direction in units of \(L_r\).
 simulation_time¶
 number_of_timesteps¶
 Duration of the simulation.
Either
simulation_time
, the simulation duration in units of \(T_r\),or
number_of_timesteps
, the total number of timesteps.
 timestep¶
 timestep_over_CFL¶
 Duration of one timestep.
Either
timestep
, in units of \(T_r\),or
timestep_over_CFL
, in units of the Courant–Friedrichs–Lewy (CFL) time.
 number_of_patches¶
A list of integers: the number of patches in each direction. Each integer must be a power of 2, and the total number of patches must be greater or equal than the number of MPI processes. It is also strongly advised to have more patches than the total number of openMP threads. See Parallelization basics.
 patch_arrangement¶
 Default
"hilbertian"
Determines the ordering of patches and the way they are separated into the various MPI processes. Options are:
"hilbertian"
: following the Hilbert curve (see this explanation)."linearized_XY"
in 2D or"linearized_XYZ"
in 3D: following the rowmajor (Cstyle) ordering."linearized_YX"
in 2D or"linearized_ZYX"
in 3D: following the columnmajor (fortranstyle) ordering. This prevents the usage of Fields diagnostics (see Parallelization basics).
 cluster_width¶
 Default
set to minimize the memory footprint of the particles pusher, especially interpolation and projection processes
For advanced users. Integer specifying the cluster width along X direction in number of cells. The “cluster” is a subpatch structure in which particles are sorted for cache improvement.
cluster_width
must divide the number of cells in one patch (in dimension X). The finest sorting is achieved withcluster_width=1
and no sorting withcluster_width
equal to the full size of a patch along dimension X. The cluster size in dimension Y and Z is always the full extent of the patch.
 maxwell_solver¶
 Default
‘Yee’
The solver for Maxwell’s equations. Only
"Yee"
and"M4"
are available for all geometries at the moment."Cowan"
,"Grassi"
,"Lehe"
and"Bouchard"
are available for2DCartesian
."Lehe"
and"Bouchard"
is available for3DCartesian
. The M4 solver is described in this paper. The Lehe solver is described in this paper. The Bouchard solver is described in this thesis p. 109
 solve_poisson¶
 Default
True
Decides if Poisson correction must be applied or not initially.
 poisson_max_iteration¶
 Default
50000
Maximum number of iteration for the Poisson solver.
 poisson_max_error¶
 Default
1e14
Maximum error for the Poisson solver.
 solve_relativistic_poisson¶
 Default
False
Decides if relativistic Poisson problem must be solved for at least one species. See Field initialization for relativistic species for more details.
 relativistic_poisson_max_iteration¶
 Default
50000
Maximum number of iteration for the Poisson solver.
 relativistic_poisson_max_error¶
 Default
1e22
Maximum error for the Poisson solver.
 EM_boundary_conditions¶
 Type
list of lists of strings
 Default
[["periodic"]]
The boundary conditions for the electromagnetic fields. Each boundary may have one of the following conditions:
"periodic"
,"silvermuller"
,"reflective"
,"ramp??"
or"PML"
.Syntax 1:[[bc_all]]
, identical for all boundaries.Syntax 2:[[bc_X], [bc_Y], ...]
, different depending on x, y or z.Syntax 3:[[bc_Xmin, bc_Xmax], ...]
, different on each boundary."silvermuller"
is an open boundary condition. The incident wave vector \(k_{inc}\) on each face is defined by"EM_boundary_conditions_k"
. When using"silvermuller"
as an injecting boundary, make sure \(k_{inc}\) is aligned with the wave you are injecting. When using"silvermuller"
as an absorbing boundary, the optimal wave absorption on a given face will be along \(k_{abs}\) the specular reflection of \(k_{inc}\) on the considered face."ramp??"
is a basic, open boundary condition designed for the spectral solver inAMcylindrical
geometry. The??
is an integer representing a number of cells (smaller than the number of ghost cells). Over the first half, the fields remain untouched. Over the second half, all fields are progressively reduced down to zero."PML"
stands for Perfectly Matched Layer. It is an open boundary condition. The number of cells in the layer is defined by"number_of_pml_cells"
. It supports laser injection as in"silvermuller"
.
Warning
In the current release, in order to use PML all
"EM_boundary_conditions"
of the simulation must be"PML"
.
 EM_boundary_conditions_k¶
 Type
list of lists of floats
 Default
[[1.,0.],[1.,0.],[0.,1.],[0.,1.]]
in 2D Default
[[1.,0.,0.],[1.,0.,0.],[0.,1.,0.],[0.,1.,0.],[0.,0.,1.],[0.,0.,1.]]
in 3D
For
silvermuller
absorbing boundaries, the x,y,z coordinates of the unit wave vectork
incident on each face (sequentially Xmin, Xmax, Ymin, Ymax, Zmin, Zmax). The number of coordinates is equal to the dimension of the simulation. The number of given vectors must be equal to 1 or to the number of faces which is twice the dimension of the simulation. In cylindrical geometry,k
coordinates are given in thexr
frame and only the Rmax face is affected.Syntax 1:[[1,0,0]]
, identical for all boundaries.Syntax 2:[[1,0,0],[1,0,0], ...]
, different on each boundary.
 number_of_pml_cells¶
 Type
list of lists of integer
 Default
[[6,6],[6,6],[6,6]]
Defines the number of cells in the
"PML"
layers using the same alternative syntaxes as"EM_boundary_conditions"
.
 time_fields_frozen¶
 Default
Time, at the beginning of the simulation, during which fields are frozen.
 reference_angular_frequency_SI¶
The value of the reference angular frequency \(\omega_r\) in SI units, only needed when collisions, ionization, radiation losses or multiphoton BreitWheeler pair creation are requested. This frequency is related to the normalization length according to \(L_r\omega_r = c\) (see Units).
 print_every¶
Number of timesteps between each info output on screen. By default, 10 outputs per simulation.
 print_expected_disk_usage¶
 Default
True
If
False
, the calculation of the expected disk usage, that is usually printed in the standard output, is skipped. This might be useful in rare cases where this calculation is costly.
 random_seed¶
 Default
0
The value of the random seed. Each patch has its own random number generator, with a seed equal to
random_seed
+ the index of the patch.
 number_of_AM¶
 Type
integer
 Default
2
The number of azimuthal modes used for the Fourier decomposition in
"AMcylindrical"
geometry. The modes range from mode 0 to mode"number_of_AM1"
.
 number_of_AM_classical_Poisson_solver¶
 Default
1
The number of azimuthal modes used for the field initialization with non relativistic Poisson solver in
"AMcylindrical"
geometry. Note that this number must be lower or equal to the number of modes of the simulation.
 number_of_AM_relativistic_field_initialization¶
 Default
1
The number of azimuthal modes used for the relativistic field initialization in
"AMcylindrical"
geometry. Note that this number must be lower or equal to the number of modes of the simulation.
 custom_oversize¶
 Type
integer
 Default
2
The number of ghostcell for each patches. The default value is set accordingly with the
interpolation_order
value.
Load Balancing¶
Load balancing (explained here) consists in exchanging
patches (domains of the simulation box) between MPI processes to reduce the
computational load imbalance.
The block LoadBalancing
is optional. If you do not define it, load balancing will
occur every 150 iterations.
LoadBalancing(
initial_balance = True,
every = 150,
cell_load = 1.,
frozen_particle_load = 0.1
)
 initial_balance¶
 Default
True
Decides if the load must be balanced at initialization. If not, the same amount of patches will be attributed to each MPI rank.
 every¶
 Default
150
Number of timesteps between each load balancing or a time selection. The value
0
suppresses all load balancing.
 cell_load¶
 Default
Computational load of a single grid cell considered by the dynamic load balancing algorithm. This load is normalized to the load of a single particle.
 frozen_particle_load¶
 Default
0.1
Computational load of a single frozen particle considered by the dynamic load balancing algorithm. This load is normalized to the load of a single particle.
Multiple decomposition of the domain¶
The block MultipleDecomposition
is necessary for spectral solvers and optional in all other cases.
When present, it activates
the Singledomain multiple decompositions (SDMD) technique
which separates the decomposition of the field grids from that of the particles.
Fields are set on large subdomain called regions (1 region per MPI process) while
particles are kept as small patches as in the standard decomposition (many patches per MPI process).
Benefits of this option are illustrated in this paper.
MultipleDecomposition(
region_ghost_cells = 2
)
 region_ghost_cells¶
 Type
integer
 Default
2
The number of ghost cells for each region. The default value is set accordingly with the
interpolation_order
. The same number of ghost cells is used in all dimensions except for spectral solver in AM geometry for which the number of radial ghost cells is always automatically set to be the same as patches.
Vectorization¶
The block Vectorization
is optional.
It controls the SIMD operations that can enhance the performance of some computations.
The technique is detailed in Ref. [Beck2019] and summarized in this doc.
It requires additional compilation options to be actived.
Vectorization(
mode = "adaptive",
reconfigure_every = 20,
initial_mode = "on"
)
 mode¶
 Default
"off"
"off"
: nonvectorized operators are used. Recommended when the number of particles per cell stays below 10."on"
: vectorized operators are used. Recommended when the number of particles per cell stays above 10. Particles are sorted per cell."adaptive"
: the best operators (scalar or vectorized) are determined and configured dynamically and locally (per patch and per species). For the moment this mode is only supported in3Dcartesian
geometry. Particles are sorted per cell.
In the
"adaptive"
mode,cluster_width
is set to the maximum.
 reconfigure_every¶
 Default
20
The number of timesteps between each dynamic reconfiguration of the vectorized operators, when using the
"adaptive"
vectorization mode. It may be set to a time selection as well.
 initial_mode¶
 Default
off
Default state when the
"adaptive"
mode is activated and no particle is present in the patch.
Moving window¶
The simulated domain can move relatively to its the initial position. The “moving window”
is (almost) periodically shifted in the x_max
direction.
Each “shift” consists in removing a column of patches from the x_min
border and
adding a new one after the x_max
border, thus changing the physical domain that the
simulation represents but keeping the same box size. This is particularly useful to
follow waves or plasma moving at high speed.
The frequency of the shifts is adjusted so that the average displacement velocity
over many shifts matches the velocity given by the user.
The user may ask for a given number of additional shifts at a given time.
These additional shifts are not taken into account for the evaluation of the average
velocity of the moving window.
The block MovingWindow
is optional. The window does not move it you do not define it.
Warning
When the window starts moving, all laser injections via SilverMuller boundary conditions are immediately stopped for physical correctness.
MovingWindow(
time_start = 0.,
velocity_x = 1.,
number_of_additional_shifts = 0.,
additional_shifts_time = 0.,
)
 time_start¶
 Type
Float.
 Default
The time at which the window starts moving.
 velocity_x¶
 Type
Float.
 Default
The average velocity of the moving window in the
x_max
direction. It muste be between 0 and 1.
 number_of_additional_shifts¶
 Type
Integer.
 Default
The number of additional shifts of the moving window.
 additional_shifts_time¶
 Type
Float.
 Default
The time at which the additional shifts are done.
Note
The particle binning diagnostics accept an “axis” called moving_x
corresponding to the x
coordinate corrected by the moving window’s current movement.
Current filtering¶
The present version of Smilei provides a multipass binomial filter on the current densities, which parameters are controlled in the following block:
CurrentFilter(
model = "binomial",
passes = [0],
kernelFIR = [0.25,0.5,0.25]
)
 model¶
 Default
"binomial"
The model for current filtering.
"binomial"
for a binomial filter."customFIR"
for a custom FIR kernel.
 passes¶
 Type
A python list of integers.
 Default
[0]
The number of passes (at each timestep) given for each dimension. If the list is of length 1, the same number of passes is assumed for all dimensions.
 kernelFIR¶
 Default
"[0.25,0.5,0.25]"
The FIR kernel for the
"customFIR"
model. The number of coefficients must be less than twice the number of ghost cells (adjusted usingcustom_oversize
).
Field filtering¶
The present version of Smilei provides a method for field filtering (at the moment, only the Friedman electric field timefilter is available) which parameters are controlled in the following block:
FieldFilter(
model = "Friedman",
theta = 0.,
)
 model¶
 Default
"Friedman"
The model for field filtering. Presently, only
"Friedman"
field filtering is available.
 theta¶
 Default
0.
The \(\theta\) parameter (between 0 and 1) of Friedman’s method.
Species¶
Each species has to be defined in a Species
block:
Species(
name = "electrons1",
position_initialization = "random",
momentum_initialization = "maxwelljuettner",
regular_number = [],
particles_per_cell = 100,
mass = 1.,
atomic_number = None,
#maximum_charge_state = None,
number_density = 10.,
# charge_density = None,
charge = 1.,
mean_velocity = [0.],
temperature = [1e10],
boundary_conditions = [
["reflective", "reflective"],
# ["periodic", "periodic"],
# ["periodic", "periodic"],
],
# thermal_boundary_temperature = None,
# thermal_boundary_velocity = None,
time_frozen = 0.0,
# ionization_model = "none",
# ionization_electrons = None,
# ionization_rate = None,
is_test = False,
pusher = "boris",
# Radiation reaction, for particles only:
radiation_model = "none",
radiation_photon_species = "photon",
radiation_photon_sampling = 1,
radiation_photon_gamma_threshold = 2,
radiation_max_emissions = 10,
# Relativistic field initialization:
relativistic_field_initialization = "False",
# For photon species only:
multiphoton_Breit_Wheeler = ["electron","positron"],
multiphoton_Breit_Wheeler_sampling = [1,1]
# Merging
merging_method = "vranic_spherical",
merge_every = 5,
merge_min_particles_per_cell = 16,
merge_max_packet_size = 4,
merge_min_packet_size = 4,
merge_momentum_cell_size = [16,16,16],
)
 name¶
The name you want to give to this species. It should be more than one character and can not start with
"m_"
.
 position_initialization¶
The method for initialization of particle positions. Options are:
"regular"
for regularly spaced. Seeregular_number
."random"
for randomly distributed."centered"
for centered in each cell (not supported inAMcylindrical
geometry.The
name
of another species from which the positions are copied. The source species must have positions initialized using one of the three other options above, and must be defined before this species.A numpy array or an HDF5 file defining all the positions of the particles. In this case you must also provide the weight of each particle (see Macroparticle weights). See Initialize particles from an array or a file.
 regular_number¶
 Type
A list of as many integers as the simulation dimension
When
position_initialization = "regular"
, this sets the number of evenlyspaced particles per cell in each direction:[Nx, Ny, Nz]
in cartesian geometries and[Nx, Nr, Ntheta]
inAMcylindrical
in which case we recommendNtheta
\(\geq 4\times (\)number_of_AM
\(1)\). If unset,particles_per_cell
must be a power of the simulation dimension, for instance, a power of 2 in2Dcartesian
.
 momentum_initialization¶
The method for initialization of particle momenta. Options are:
"maxwelljuettner"
for a relativistic maxwellian (see how it is done)"rectangular"
for a rectangular distribution"cold"
for zero temperatureA numpy array or an HDF5 file defining all the momenta of the particles. See Initialize particles from an array or a file.
The first 2 distributions depend on the parameter
temperature
explained below.
 mass¶
The mass of particles, in units of the electron mass \(m_e\).
 atomic_number¶
 Default
0
The atomic number of the particles, required only for ionization. It must be lower than 101.
 maximum_charge_state¶
 Default
0
The maximum charge state of a species for which the ionization model is
"from_rate"
.
 number_density¶
 charge_density¶
 Type
float or profile
The absolute value of the charge density or number density (choose one only) of the particle distribution, in units of the reference density \(N_r\) (see Units).
 mean_velocity¶
 Type
a list of 3 floats or profiles
The initial drift velocity of the particles, in units of the speed of light \(c\).
WARNING: For massless particles, this is actually the momentum in units of \(m_e c\).
 temperature¶
 Type
a list of 3 floats or profiles
The initial temperature of the particles, in units of \(m_ec^2\).
 boundary_conditions¶
 Type
a list of lists of strings
 Default
[["periodic"]]
The boundary conditions for the particles of this species. Each boundary may have one of the following conditions:
"periodic"
,"reflective"
,"remove"
(particles are deleted),"stop"
(particle momenta are set to 0), and"thermalize"
. For photon species (mass=0
), the last two options are not available.Syntax 1:[[bc_all]]
, identical for all boundaries.Syntax 2:[[bc_X], [bc_Y], ...]
, different depending on x, y or z.Syntax 3:[[bc_Xmin, bc_Xmax], ...]
, different on each boundary.
 thermal_boundary_temperature¶
 Default
None
A list of floats representing the temperature of the thermal boundaries (those set to
"thermalize"
inboundary_conditions
) for each spatial coordinate. Currently, only the first coordinate (x) is taken into account.
 thermal_boundary_velocity¶
 Default
[]
A list of floats representing the components of the particles’ drift velocity after encountering the thermal boundaries (those set to
"thermalize"
inboundary_conditions
).
 time_frozen¶
 Default
The time during which the particles are “frozen”, in units of \(T_r\). Frozen particles do not move and therefore do not deposit any current density either. Nonetheless, they deposit a charge density. They are computationally much cheaper than nonfrozen particles and oblivious to any EMfields in the simulation. Note that frozen particles can be ionized (this is computationally much cheaper if ion motion is not relevant).
 ionization_model¶
 Default
"none"
The model for ionization:
"tunnel"
for field ionization (requires species with anatomic_number
)"tunnel_envelope_averaged"
for field ionization with a laser envelope"from_rate"
, relying on a userdefined ionization rate (requires species with amaximum_charge_state
).
 ionization_rate¶
A python function giving the userdefined ionisation rate as a function of various particle attributes. To use this option, the numpy package must be available in your python installation. The function must have one argument, that you may call, for instance,
particles
. This object has several attributesx
,y
,z
,px
,py
,pz
,charge
,weight
andid
. Each of these attributes are provided as numpy arrays where each cell corresponds to one particle.The following example defines, for a species with maximum charge state of 2, an ionization rate that depends on the initial particle charge and linear in the x coordinate:
from numpy import exp, zeros_like def my_rate(particles): rate = zeros_like(particles.x) charge_0 = (particles.charge==0) charge_1 = (particles.charge==1) rate[charge_0] = r0 * particles.x[charge_0] rate[charge_1] = r1 * particles.x[charge_1] return rate Species( ..., ionization_rate = my_rate )
 ionization_electrons¶
The name of the electron species that
ionization_model
uses when creating new electrons.
 is_test¶
 Default
False
Flag for test particles. If
True
, this species will contain only test particles which do not participate in the charge and currents.
 pusher¶
 Default
"boris"
Type of pusher to be used for this species. Options are:
"boris"
: The relativistic Boris pusher"borisnr"
: The nonrelativistic Boris pusher"vay"
: The relativistic pusher of J. L. Vay"higueracary"
: The relativistic pusher of A. V. Higuera and J. R. Cary"norm"
: For photon species only (rectilinear propagation)"ponderomotive_boris"
: modified relativistic Boris pusher for species interacting with the laser envelope model. Valid only if the species has nonzero mass
 radiation_model¶
 Default
"none"
The radiation reaction model used for this species (see Highenergy photon emission & radiation reaction).
"none"
: no radiation"LandauLifshitz"
(orll
): LandauLifshitz model approximated for high energies"correctedLandauLifshitz"
(orcll
): with quantum correction"Niel"
: a stochastic radiation model based on the work of Niel et al.."MonteCarlo"
(ormc
): MonteCarlo radiation model. This model can be configured to generate macrophotons withradiation_photon_species
.
This parameter cannot be assigned to photons (mass = 0).
Radiation is emitted only with the
"MonteCarlo"
model whenradiation_photon_species
is defined.
 radiation_photon_species¶
The
name
of the photon species in which the MonteCarloradiation_model
will generate macrophotons. If unset (orNone
), no macrophoton will be created. The target photon species must be have its mass set to 0, and appear after the particle species in the namelist.This parameter cannot be assigned to photons (mass = 0).
 radiation_photon_sampling¶
 Default
1
The number of macrophotons generated per emission event, when the macrophoton creation is activated (see
radiation_photon_species
). The total macrophoton weight is still conserved.A large number may rapidly slow down the performances and lead to memory saturation.
This parameter cannot be assigned to photons (mass = 0).
 radiation_max_emissions¶
 Default
10
The maximum number of emission MonteCarlo event a macroparticle can undergo during a timestep. Since this value is used to allocate some buffers, a high value can saturate memory.
This parameter cannot be assigned to photons (mass = 0).
 radiation_photon_gamma_threshold¶
 Default
2
The threshold on the photon energy for the macrophoton emission when using the radiation reaction MonteCarlo process. Under this threshold, the macrophoton from the radiation reaction MonteCarlo process is not created but still taken into account in the energy balance. The default value corresponds to twice the electron rest mass energy that is the required energy to decay into electronpositron pairs.
This parameter cannot be assigned to photons (mass = 0).
 relativistic_field_initialization¶
 Default
False
Flag for relativistic particles. If
True
, the electromagnetic fields of this species will added to the electromagnetic fields already present in the simulation. This operation will be performed when time equalstime_frozen
. See Field initialization for relativistic species for details on the computation of the electromagentic fields of a relativistic species. To have physically meaningful results, we recommend to place a species which requires this method of field initialization far from other species, otherwise the latter could experience instantly turnedon unphysical forces by the relativistic species’ fields.
 multiphoton_Breit_Wheeler¶
 Default
[None,None]
An list of the
name
of two species: electrons and positrons created through the Multiphoton BreitWheeler pair creation. By default, the process is not activated.This parameter can only be assigned to photons species (mass = 0).
 multiphoton_Breit_Wheeler_sampling¶
 Default
[1,1]
A list of two integers: the number of electrons and positrons generated per photon decay in the Multiphoton BreitWheeler pair creation. The total macroparticle weight is still conserved.
Large numbers may rapidly slow down the performances and lead to memory saturation.
This parameter can only be assigned to photons species (mass = 0).
Particle Injector¶
Injectors enable to inject macroparticles in the simulation domain from the boundaries.
By default, some parameters that are not specified are inherited from the associated species
.
Each particle injector has to be defined in a ParticleInjector
block:
ParticleInjector(
name = "injector1",
species = "electrons1",
box_side = "xmin",
time_envelope = tgaussian(start=0, duration=10., order=4),
# Parameters inherited from the associated ``species`` by default
position_initialization = "species",
momentum_initialization = "rectangular",
mean_velocity = [0.5,0.,0.],
temperature = [1e30],
number_density = 1,
particles_per_cell = 16,
)
 name¶
The name you want to give to this injector. If you do not specify a name, it will be attributed automatically. The name is useful if you want to inject particles at the same position of another injector.
 species¶
The name of the species in which to inject the new particles
 box_side¶
From where the macroparticles are injected. Options are:
"xmin"
"xmax"
"ymin"
"ymax"
"zmax"
"zmin"
 time_envelope¶
 Type
a python function or a time profile
 Default
tconstant()
The temporal envelope of the injector.
 position_initialization¶
 Default
parameters provided the species
The method for initialization of particle positions. Options are:
"species"
or empty""
: injector uses the option of the specifiedspecies
."regular"
for regularly spaced. Seeregular_number
."random"
for randomly distributed"centered"
for centered in each cellThe
name
of another injector from which the positions are copied. This option requires (1) that the target injector’s positions are initialized using one of the three other options above.
 momentum_initialization¶
 Default
parameters provided the species
The method for initialization of particle momenta. Options are:
"species"
or empty""
: injector uses the option of the specifiedspecies
."maxwelljuettner"
for a relativistic maxwellian (see how it is done)"rectangular"
for a rectangular distribution
 mean_velocity¶
 Type
a list of 3 floats or profiles
 Default
parameters provided the species
The initial drift velocity of the particles, in units of the speed of light \(c\).
WARNING: For massless particles, this is actually the momentum in units of \(m_e c\).
 temperature¶
 Type
a list of 3 floats or profiles
 Default
parameters provided the species
The initial temperature of the particles, in units of \(m_ec^2\).
 particles_per_cell¶
 Type
float or profile
 Default
parameters provided the species
The number of particles per cell to use for the injector.
 number_density¶
 charge_density¶
 Type
float or profile
 Default
parameters provided the species
The absolute value of the number density or charge density (choose one only) of the particle distribution, in units of the reference density \(N_r\) (see Units)
 regular_number¶
 Type
A list of as many integers as the simulation dimension
Same as for Species. When
position_initialization = "regular"
, this sets the number of evenlyspaced particles per cell in each direction:[Nx, Ny, Nz]
in cartesian geometries.
Particle Merging¶
The macroparticle merging method is documented in
the corresponding page.
Note that for merging to be able to operate either vectorization or cell sorting must be activated.
It is optionnally specified in the Species
block:
Species(
....
# Merging
merging_method = "vranic_spherical",
merge_every = 5,
merge_min_particles_per_cell = 16,
merge_max_packet_size = 4,
merge_min_packet_size = 4,
merge_momentum_cell_size = [16,16,16],
merge_discretization_scale = "linear",
# Extra parameters for experts:
merge_min_momentum_cell_length = [1e10, 1e10, 1e10],
merge_accumulation_correction = True,
)
 merging_method¶
 Default
"none"
The particle merging method to use:
"none"
: no merging"vranic_cartesian"
: method of M. Vranic with a cartesian momentumspace decomposition"vranic_spherical"
: method of M. Vranic with a spherical momentumspace decomposition
 merge_every¶
 Default
0
Number of timesteps between each merging event or a time selection.
 min_particles_per_cell¶
 Default
4
The minimum number of particles per cell for the merging.
 merge_min_packet_size¶
 Default
4
The minimum number of particles per packet to merge. Must be greater or equal to 4.
 merge_max_packet_size¶
 Default
4
The maximum number of particles per packet to merge.
 merge_momentum_cell_size¶
 Default
[16,16,16]
A list of 3 integers defining the number of subgroups in each direction for the momentumspace discretization.
 merge_discretization_scale¶
 Default
"linear"
The momentum discretization scale::
"linear"
or"log"
. The"log"
scale only works with the spherical discretization at the moment.
 merge_min_momentum¶
 Default
1e5
[for experts] The minimum momentum value when the log scale is chosen (
merge_discretization_scale = log
). This avoids a potential 0 value in the log domain.
 merge_min_momentum_cell_length¶
 Default
[1e10,1e10,1e10]
[for experts] The minimum subgroup length for the momentumspace discretization (below which the number of subgroups is set to 1).
 merge_accumulation_correction¶
 Default
True
[for experts] Activates the accumulation correction (see Particle Merging for more information). The correction only works in linear scale.
Lasers¶
A laser consists in applying oscillating boundary conditions for the magnetic
field on one of the box sides. The only boundary condition that supports lasers
is "silvermuller"
(see EM_boundary_conditions
).
There are several syntaxes to introduce a laser in Smilei:
Note
The following definitions are given for lasers incoming from the xmin
or xmax
boundaries. For lasers incoming from ymin
or ymax
, replace the By
profiles by Bx
profiles. For lasers incoming from zmin
or zmax
,
replace By
and Bz
profiles by Bx
and By
profiles, respectively.
1. Defining a generic wave
Laser( box_side = "xmin", space_time_profile = [ By_profile, Bz_profile ] space_time_profile_AM = [ Br_mode0, Bt_mode0, Br_mode1, Bt_mode1, ... ] )
 box_side¶
 Default
"xmin"
Side of the box from which the laser originates:
"xmin"
,"xmax"
,"ymin"
,"ymax"
,"zmin"
or"zmax"
.In the cases of
"ymin"
or"ymax"
, replace, in the following profiles, coordinates y by x, and fields \(B_y\) by \(B_x\).In the cases of
"zmin"
or"zmax"
, replace, in the following profiles, coordinates y by x, coordinates z by y, fields \(B_y\) by \(B_x\) and fields \(B_z\) by \(B_y\).
 space_time_profile¶
 Type
A list of two python functions
The full wave expression at the chosen box side. It is a list of two python functions taking several arguments depending on the simulation dimension: \((t)\) for a 1D simulation, \((y,t)\) for a 2D simulation (etc.) The two functions represent \(B_y\) and \(B_z\), respectively. This can be used only in Cartesian geometries.
 space_time_profile_AM¶
 Type
A list of maximum 2 x
number_of_AM
python functions.
These profiles define the first modes of \(B_r\) and \(B_\theta\) in the order shown in the above example. Undefined modes are considered zero. This can be used only in
AMcylindrical
geometry. In this geometry a twodimensional \((x,r)\) grid is used and the laser is injected from a \(x\) boundary, thus the provided profiles must be a function of \((r,t)\).
2. Defining the wave envelopes
Laser( box_side = "xmin", omega = 1., chirp_profile = tconstant(), time_envelope = tgaussian(), space_envelope = [ By_profile , Bz_profile ], phase = [ PhiY_profile, PhiZ_profile ], delay_phase = [ 0., 0. ] )This implements a wave of the form:
\[ \begin{align}\begin{aligned}B_y(\mathbf{x}, t) = S_y(\mathbf{x})\; T\left(tt_{0y}\right) \;\sin\left( \omega(t) t  \phi_y(\mathbf{x}) \right)\\B_z(\mathbf{x}, t) = S_z(\mathbf{x})\; T\left(tt_{0z}\right) \;\sin\left( \omega(t) t  \phi_z(\mathbf{x}) \right)\end{aligned}\end{align} \]where \(T\) is the temporal envelope, \(S_y\) and \(S_z\) are the spatial envelopes, \(\omega\) is the timevarying frequency, \(\phi_y\) and \(\phi_z\) are the phases, and we defined the delays \(t_{0y} = (\phi_y(\mathbf{x})\varphi_y)/\omega(t)\) and \(t_{0z} = (\phi_z(\mathbf{x})\varphi_z)/\omega(t)\).
 omega¶
 Default
The laser angular frequency.
 chirp_profile¶
 Type
a python function or a time profile
 Default
tconstant()
The variation of the laser frequency over time, such that \(\omega(t)=\)
omega
xchirp_profile
\((t)\).Warning
This definition of the chirp profile is not standard. Indeed, \(\omega(t)\) as defined here is not the instantaneous frequency, \(\omega_{\rm inst}(t)\), which is obtained from the time derivative of the phase \(\omega(t) t\).
Should one define the chirp as \(C(t) = \omega_{\rm inst}(t)/\omega\) (with \(\omega\) defined by the input parameter \(\mathtt{omega}\)), the user can easily obtain the corresponding chirp profile as defined in Smilei as:
\[\mathtt{chirp\_profile}(t) = \frac{1}{t} \int_0^t dt' C(t')\,.\]Let us give as an example the case of a linear chirp, with the instantaneous frequency \(\omega_{\rm inst}(t) = \omega [1+\alpha\,\omega(tt_0)]\). \(C(t) = 1+\alpha\,\omega(tt_0)\). The corresponding input chirp profile reads:
\[\mathtt{chirp\_profile}(t) = 1  \alpha\, \omega t_0 + \frac{\alpha}{2} \omega t\]Similarly, for a geometric (exponential) chirp such that \(\omega_{\rm inst}(t) = \omega\, \alpha^{\omega t}\), \(C(t) = \alpha^{\omega t}\), and the corresponding input chirp profile reads:
\[\mathtt{chirp\_profile}(t) = \frac{\alpha^{\omega t}  1}{\omega t \, \ln \alpha}\,.\]
 time_envelope¶
 Type
a python function or a time profile
 Default
tconstant()
The temporal envelope of the laser.
 space_envelope¶
 Type
a list of two python functions or two spatial profiles
 Default
[ 1., 0. ]
The two spatial envelopes \(S_y\) and \(S_z\).
 phase¶
 Type
a list of two python functions or two spatial profiles
 Default
[ 0., 0. ]
The two spatiallyvarying phases \(\phi_y\) and \(\phi_z\).
3. Defining a 1D planar wave
For onedimensional simulations, you may use the simplified laser creator:
LaserPlanar1D( box_side = "xmin", a0 = 1., omega = 1., polarization_phi = 0., ellipticity = 0., time_envelope = tconstant() )
 a0¶
 Default
The normalized vector potential
 polarization_phi¶
 Default
The angle of the polarization ellipse major axis relative to the XY plane, in radians.
 ellipticity¶
 Default
The polarization ellipticity: 0 for linear and \(\pm 1\) for circular.
4. Defining a 2D gaussian wave
For twodimensional simulations, you may use the simplified laser creator:
LaserGaussian2D( box_side = "xmin", a0 = 1., omega = 1., focus = [50., 40.], waist = 3., incidence_angle = 0., polarization_phi = 0., ellipticity = 0., time_envelope = tconstant() )
 focus¶
 Type
A list of two floats
[X, Y]
The
X
andY
positions of the laser focus.
 waist¶
The waist value. Transverse coordinate at which the field is at 1/e of its maximum value.
 incidence_angle¶
 Default
The angle of the laser beam relative to the normal to the injection plane, in radians.
 time_envelope¶
Time envelope of the field (not intensity).
5. Defining a 3D gaussian wave
For threedimensional simulations, you may use the simplified laser creator:
LaserGaussian3D( box_side = "xmin", a0 = 1., omega = 1., focus = [50., 40., 40.], waist = 3., incidence_angle = [0., 0.1], polarization_phi = 0., ellipticity = 0., time_envelope = tconstant() )This is almost the same as
LaserGaussian2D
, with thefocus
parameter having now 3 elements (focus position in 3D), and theincidence_angle
being a list of two angles, corresponding to rotations aroundy
andz
, respectively.When injecting on
"ymin"
or"ymax"
, the incidence angles corresponds to rotations aroundx
andz
, respectively.
6. Defining a gaussian wave with Azimuthal Fourier decomposition
For simulations with
"AMcylindrical"
geometry, you may use the simplified laser creator:LaserGaussianAM( box_side = "xmin", a0 = 1., omega = 1., focus = [50., 0.], waist = 3., polarization_phi = 0., ellipticity = 0., time_envelope = tconstant() )Note that here, the focus is given in [x,r] coordinates.
7. Defining a generic wave at some distance from the boundary
In some cases, the laser field is not known at the box boundary, but rather at some plane inside the box. Smilei can precalculate the corresponding wave at the boundary using the angular spectrum method. This technique is only available in 2D and 3D cartesian geometries and requires the python packages numpy. A detailed explanation of the method is available. The laser is introduced using:
LaserOffset( box_side = "xmin", space_time_profile = [ By_profile, Bz_profile ], offset = 10., extra_envelope = tconstant(), keep_n_strongest_modes = 100, angle = 10./180.*3.14159 )
 space_time_profile¶
 Type
A list of two python functions
The magnetic field profiles at some arbitrary plane, as a function of space and time. The arguments of these profiles are
(y,t)
in 2D and(y,z,t)
in 3D.
 offset¶
The distance from the box boundary to the plane where
space_time_profile
is defined.
 extra_envelope¶
 Type
a python function or a python profile
 Default
lambda *z: 1.
, which means a profile of value 1 everywhereAn extra envelope applied at the boundary, on top of the
space_time_profile
. This envelope takes two arguments (y
,t
) in 2D, and three arguments (y
,z
,t
) in 3D. As the wave propagation technique stores a limited number of Fourier modes (in the time domain) of the wave, some periodicity can be obtained in the actual laser. One may thus observe that the laser pulse is repeated several times. The envelope can be used to remove these spurious repetitions.
 keep_n_strongest_modes¶
 Default
100
The number of temporal Fourier modes that are kept during the preprocessing. See this page for more details.
 angle¶
 Default
Angle between the boundary and the profile’s plane, the rotation being around \(z\). See this page for more details.
 fft_time_window¶
 Default
Time during which the
space_time_profile
is sampled (calculating theLaserOffset
on the whole simulation duration can be costly). Note that the Fourier approach will naturally repeat the signal periodically.
 fft_time_step¶
 Default
Temporal step between each sample of the
space_time_profile
. Chosing a larger step can help reduce the memory load but will remove high temporal frequencies.
 number_of_processes¶
 Default
all available processes
The number of MPI processes that will be used for computing the
LaserOffset
. Using more processes computes the FFT faster, but too many processes may be very costly in communication. In addition, using too few may not allow the arrays to fit in memory.
 file¶
 Default
None
The path to a
LaserOffset*.h5
file generated from a previous simulation. This option can help reduce the computation time by reusing theLaserOffset
computation from a previous simulation.
Laser envelope model¶
In all the available geometries, it is possible to model a laser pulse
propagating in the x
direction
using an envelope model (see Laser envelope model for the advantages
and limits of this approximation).
The fast oscillations of the laser are neglected and all the physical
quantities of the simulation, including the electromagnetic fields and
their source terms, as well as the particles positions and momenta, are
meant as an average over one or more optical cycles.
Effects involving characteristic lengths comparable to the laser central
wavelength (i.e. sharp plasma density profiles) cannot be modeled with
this option.
Note
The envelope model in "AMcylindrical"
geometry is implemented only in the hypothesis of
cylindrical symmetry, i.e. only one azimuthal mode. Therefore, to use it the user must choose
number_of_AM = 1
.
Contrarily to a standard Laser initialized with the SilverMüller boundary conditions, the laser envelope will be entirely initialized inside the simulation box at the start of the simulation.
Currently only one laser pulse of a given frequency propagating in the positive x direction can be speficified. However, a multipulse setup can be initialized if a multipulse profile is specified, e.g. if the temporal profile is given by two adjacents gaussian functions. The whole multipulse profile would have the same carrier frequency and would propagate in the positive x direction. For the moment it is not possible to specify more than one laser envelope profile, e.g. two counterpropagating lasers, or two lasers with different carrier frequency.
Please note that describing a laser through its complex envelope loses physical accuracy if its characteristic spacetime variation scales are too small, i.e. of the order of the laser central wavelength (see Laser envelope model). Thus, spacetime profiles with variation scales larger than this length should be used.
1. Defining a generic laser envelope
Following is the generic laser envelope creator
LaserEnvelope(
omega = 1.,
envelope_solver = 'explicit',
envelope_profile = envelope_profile,
Envelope_boundary_conditions = [["reflective"]]
polarization_phi = 0.,
ellipticity = 0.
)
 omega¶
 Default
1.
The laser angular frequency.
 envelope_profile¶
 Type
a python function or a python profile
 Default
None
The laser spacetime profile, so if the geometry is
3Dcartesian
a function of 4 arguments (3 for space, 1 for time) is necessary. Please note that the envelope will be entirely initialized in the simulation box already at the start of the simulation, so the time coordinate will be applied to thex
direction instead of time. It is recommended to initialize the laser envelope in vacuum, separated from the plasma, to avoid unphysical results. Envelopes with variation scales near to the laser wavelength do not satisfy the assumptions of the envelope model (see Laser envelope model), yielding inaccurate results.
 envelope_solver¶
 Default
explicit
The solver scheme for the envelope equation.
"explicit"
: an explicit scheme based on central finite differences."explicit_reduced_dispersion"
: the finite difference derivatives alongx
in the"explicit"
solver are substituted by optimized derivatives to reduce numerical dispersion.
 Envelope_boundary_conditions¶
 Type
list of lists of strings
 Default
[["reflective"]]
For the moment, only reflective boundary conditions are implemented in the resolution of the envelope equation.
 polarization_phi¶
 Default
The angle of the polarization ellipse major axis relative to the XY plane, in radians. Needed only for ionization.
 ellipticity¶
 Default
The polarization ellipticity: 0 for linear and 1 for circular. For the moment, only these two polarizations are available.
2. Defining a 1D laser envelope
Following is the simplified laser envelope creator in 1D
LaserEnvelopePlanar1D(
a0 = 1.,
time_envelope = tgaussian(center=150., fwhm=40.),
envelope_solver = 'explicit',
Envelope_boundary_conditions = [ ["reflective"] ],
polarization_phi = 0.,
ellipticity = 0.
)
3. Defining a 2D gaussian laser envelope
Following is the simplified gaussian laser envelope creator in 2D
LaserEnvelopeGaussian2D(
a0 = 1.,
focus = [150., 40.],
waist = 30.,
time_envelope = tgaussian(center=150., fwhm=40.),
envelope_solver = 'explicit',
Envelope_boundary_conditions = [ ["reflective"] ],
polarization_phi = 0.,
ellipticity = 0.
)
4. Defining a 3D gaussian laser envelope
Following is the simplified laser envelope creator in 3D
LaserEnvelopeGaussian3D(
a0 = 1.,
focus = [150., 40., 40.],
waist = 30.,
time_envelope = tgaussian(center=150., fwhm=40.),
envelope_solver = 'explicit',
Envelope_boundary_conditions = [ ["reflective"] ],
polarization_phi = 0.,
ellipticity = 0.
)
5. Defining a cylindrical gaussian laser envelope
Following is the simplified laser envelope creator in "AMcylindrical"
geometry (remember that
in this geometry the envelope model can be used only if number_of_AM = 1
)
LaserEnvelopeGaussianAM(
a0 = 1.,
focus = [150., 40.],
waist = 30.,
time_envelope = tgaussian(center=150., fwhm=40.),
envelope_solver = 'explicit',
Envelope_boundary_conditions = [ ["reflective"] ],
polarization_phi = 0.,
ellipticity = 0.
)
The arguments appearing LaserEnvelopePlanar1D
, LaserEnvelopeGaussian2D
,
LaserEnvelopeGaussian3D
and LaserEnvelopeGaussianAM
have the same meaning they would have in a
normal LaserPlanar1D
, LaserGaussian2D
, LaserGaussian3D
and LaserGaussianAM
,
with some differences:
 time_envelope¶
Since the envelope will be entirely initialized in the simulation box already at the start of the simulation, the time envelope will be applied in the
x
direction instead of time. It is recommended to initialize the laser envelope in vacuum, separated from the plasma, to avoid unphysical results. Temporal envelopes with variation scales near to the laser wavelength do not satisfy the assumptions of the envelope model (see Laser envelope model), yielding inaccurate results.
 waist¶
Please note that a waist size comparable to the laser wavelength does not satisfy the assumptions of the envelope model.
It is important to remember that the profile defined through the blocks
LaserEnvelopePlanar1D
, LaserEnvelopeGaussian2D
, LaserEnvelopeGaussian3D
correspond to the complex envelope of the laser vector potential component
\(\tilde{A}\) in the polarization direction.
The calculation of the correspondent complex envelope for the laser electric field
component in that direction is described in Laser envelope model.
Note that only order 2 interpolation and projection are supported in presence of the envelope model for the laser.
The parameters polarization_phi
and ellipticity
specify the polarization state of the laser. In envelope model implemented in Smilei,
they are only used to compute the rate of ionization and the initial momentum of the electrons newly created by ionization,
where the polarization of the laser plays an important role (see Ionization).
For all other purposes (e.g. the particles equations of motions, the computation of the ponderomotive force,
the evolution of the laser), the polarization of the laser plays no role in the envelope model.
External fields¶
An constant external field can be applied over the whole box
(at the beginning of the simulation) using an ExternalField
block:
ExternalField(
field = "Ex",
profile = constant(0.01, xvacuum=0.1)
)
 field¶
Field name:
"Ex"
,"Ey"
,"Ez"
,"Bx"
,"By"
or"Bz"
.
 profile¶
 Type
float or profile
The initial spatial profile of the applied field. Refer to Units to understand the units of this field.
Prescribed fields¶
Userdefined electromagnetic fields, with spatiotemporal dependence, can be superimposed to the selfconsistent Maxwell fields. These fields push the particles but do not participate in the Maxwell solver: they are not selfconsistent. They are however useful to describe charged particles’ dynamics in a given electromagnetic field.
This feature is accessible using the PrescribedField
block:
from numpy import cos, sin
def myPrescribedProfile(x,t):
return cos(x)*sin(t)
PrescribedField(
field = "Ex",
profile = myPrescribedProfile
)
 field¶
Field name:
"Ex"
,"Ey"
,"Ez"
,"Bx_m"
,"By_m"
or"Bz_m"
.
Warning
When prescribing a magnetic field, always use the timecentered fields "Bx_m"
, "By_m"
or "Bz_m"
.
These fields are those used in the particle pusher, and are defined at integer timesteps.
 profile¶
 Type
float or profile
The spatiotemporal profile of the applied field: a python function with arguments (x, t) or (x, y, t), etc. Refer to Units to understand the units of this field.
Antennas¶
An antenna is an extra current applied during the whole simulation.
It is applied using an Antenna
block:
Antenna(
field = "Jz",
space_profile = gaussian(0.01),
time_profile = tcosine(base=0., duration=1., freq=0.1)
)
 field¶
The name of the current:
"Jx"
,"Jy"
or"Jz"
.
 space_profile¶
 Type
float or profile
The initial spatial profile of the applied antenna. Refer to Units to understand the units of this current.
 time_profile¶
 Type
float or profile
The temporal profile of the applied antenna. It multiplies
space_profile
.
 space_time_profile¶
 Type
float or profile
A space & time profile for the antenna (not compatible with
space_profile
ortime_profile
). It should haveN+1``arguments, where ``N
is the dimension of the simulation. For instance(x,t)
in 1D,(x,y,t)
in 2D, etc.The function must accept
x
,y
andz
either as floats or numpy arrays. If it accepts floats, the return value must be a float. If it accepts numpy arrays, these arrays will correspond to the coordinates of 1 patch, and the return value must be a numpy array of the same size.
Walls¶
A wall can be introduced using a PartWall
block in order to
reflect, stop, thermalize or kill particles which reach it:
PartWall(
kind = "reflective",
x = 20.
)
 kind¶
The kind of wall:
"reflective"
,"stop"
,"thermalize"
or"remove"
.
Collisions & reactions¶
Binary collisions & reactions account for shortrange Coulomb interactions of particles (shorter than the cell size), but also include other effects such as impact ionization and nuclear reactions. These are gathered under this section because they are treated as binary processes (meaning they happen during the encounter of two macroparticles).
They are specified by one or several Collisions
blocks:
Collisions(
species1 = ["electrons1", "electrons2"],
species2 = ["ions1"],
debug_every = 1000,
coulomb_log = 0.,
coulomb_log_factor = 1.,
ionizing = False,
# nuclear_reaction = [],
)
 species1¶
 species2¶
Lists of species’
name
.The collisions and reactions will occur between all species under the group
species1
and all species under the groupspecies2
. For example, to collide all electrons with ions:species1 = ["electrons1", "electrons2"], species2 = ["ions"]
Warning
This does not make
electrons1
collide withelectrons2
.The two groups of species have to be completely different OR exactly equal. In other words, if
species1
is not equal tospecies2
, then they cannot have any common species. If the two groups are exactly equal, we call this situation intracollisions.Note
If both lists
species1
andspecies2
contain only one species, the algorithm is potentially faster than the situation with several species in one or the other list. This is especially true if the machine accepts SIMD vectorization.
 every¶
 Default
1
Number of timesteps between each computation of the collisions or reactions. Use a number higher than 1 only if you know the collision frequency is low with respect to the inverse of the timestep.
 debug_every¶
 Default
0
Number of timesteps between each output of information about collisions or reactions. If 0, there will be no outputs.
 time_frozen¶
 Default
The time during which no collisions or reactions happen, in units of \(T_r\).
 coulomb_log¶
 Default
The Coulomb logarithm.
If \(= 0\), the Coulomb logarithm is automatically computed for each collision.
If \(> 0\), the Coulomb logarithm is equal to this value.
If \(< 0\), collisions are not treated (but other reactions may happen).
 coulomb_log_factor¶
 Default
A constant, strictly positive factor that multiplies the Coulomb logarithm, regardless of
coulomb_log
being automatically computed or set to a constant value. This can help, for example, to compensate artificiallyreduced ion masses.
 ionizing¶
 Default
False
Collisional ionization is set when this parameter is not
False
. It can either be set to the name of a preexisting electron species (where the ionized electrons are created), or toTrue
(the first electron species inspecies1
orspecies2
is then chosen for ionized electrons).One of the species groups must be all electrons (
mass
= 1), and the other one all ions of the sameatomic_number
.
 nuclear_reaction¶
 Type
a list of strings
 Default
None
(no nuclear reaction)
A list of the species names for the products of Nuclear reactions that may occur during collisions. You may omit product species if they are not necessary for the simulation.
All members of
species1
must be the same type of atoms, which is automatically recognized by theirmass
andatomic_number
. The same applies for all members ofspecies2
.In the current version, only the reaction D(d,n)He³ is available.
 nuclear_reaction_multiplier¶
 Type
a float
 Default
(automatically adjusted)
The rate multiplier for nuclear reactions. It is a positive number that artificially increases the occurence of reactions so that a good statistics is obtained. The number of actual reaction products is adjusted by changing their weights in order to provide a physically correct number of reactions. Leave this number to
0.
for an automatic rate multiplier: the final number of produced macroparticles will be of the same order as that of reactants.
Radiation reaction¶
The block RadiationReaction()
enables to tune the radiation loss properties
(see Highenergy photon emission & radiation reaction).
Many parameters are used for the generation of the crosssection tables
for the MonteCarlo emission process.
If the tables already exist in the simulation directory, then they will be read
and no new table will be generated by Smilei.
Otherwise, Smilei can compute and output these
tables.
RadiationReaction(
# Radiation parameters
minimum_chi_continuous = 1e3,
minimum_chi_discontinuous = 1e2,
table_path = "<path to the external table folder>",
# Parameters for Niel et al.
Niel_computation_method = "table",
)
 minimum_chi_continuous¶
 Default
1e3
Threshold on the particle quantum parameter particle_chi. When a particle has a quantum parameter below this threshold, radiation reaction is not taken into account.
 minimum_chi_discontinuous¶
 Default
1e2
Threshold on the particle quantum parameter particle_chi between the continuous and the discontinuous radiation model.
 table_path¶
 Default
""
Path to the directory that contains external tables for the radiation losses. If empty, the default tables are used. Default tables are embedded in the code. External tables can be generated using the external tool smilei_tables (see Generation of the external tables).
 Niel_computation_method¶
 Default
"table"
Method to compute the value of the table h of Niel et al during the emission process. The possible values are:
"table"
: the h function is tabulated. The table is computed at initialization or read from an external file."fit5"
: A polynomial fit of order 5 is used. No table is required. The maximal relative error to the reference data is of maximum of 0.02. The fit is valid for quantum parameters \(\chi\) between 1e3 and 10."fit10"
: A polynomial fit of order 10 is used. No table is required. The precision if better than the fit of order 5 with a maximal relative error of 0.0002. The fit is valid for quantum parameters \(\chi\) between 1e3 and 10."ridgers"
: The fit of Ridgers given in Ridgers et al., ArXiv 1708.04511 (2017)
The use of tabulated values is best for accuracy but not for performance. Table access prevent total vectorization. Fits are vectorizable.
Multiphoton BreitWheeler¶
The block MultiphotonBreitWheeler
enables to tune parameters of the
multiphoton BreitWheeler process and particularly the table generation.
For more information on this physical mechanism, see Multiphoton BreitWheeler pair creation.
There are three tables used for the multiphoton BreitWheeler refers to as the integration_dT_dchi, min_particle_chi_for_xi and xi table.
MultiphotonBreitWheeler(
# Path to the tables
table_path = "<path to the external table folder>",
)
 table_path¶
 Default
""
Path to the directory that contains external tables for the multiphoton BreitWheeler. If empty, the default tables are used. Default tables are embedded in the code. External tables can be generated using the external tool smilei_tables (see Generation of the external tables).
Scalar diagnostics¶
Smilei can collect various scalar data, such as total particle energy, total field energy, etc.
This is done by including the block DiagScalar
:
DiagScalar(
every = 10 ,
vars = ["Utot", "Ukin", "Uelm"],
precision = 10
)
 every¶
Number of timesteps between each output or a time selection.
 vars¶
 Default
[]
List of scalars that will be actually output. Note that most scalars are computed anyways.Omit this argument to include all scalars.
 precision¶
 Default
10
Number of digits of the outputs.
Warning
Scalars diagnostics min/max cell are not yet supported in "AMcylindrical"
geometry.
The full list of available scalars is given in the table below.
Warning
As some of these quantities are integrated in space and/or time, their units are unusual, and depend on the simulation dimension. All details here.
Spaceintegrated energy densities 



Space & timeintegrated Energies lost/gained at boundaries 



Particle information 



Fields information 


Checkout the postprocessing documentation as well.
Fields diagnostics¶
Smilei can collect various field data (electromagnetic fields, currents and density)
taken at the location of the PIC grid, both as instantaneous values and averaged values.
This is done by including a block DiagFields
:
DiagFields(
#name = "my field diag",
every = 10,
time_average = 2,
fields = ["Ex", "Ey", "Ez"],
#subgrid = None
)
 name¶
Optional name of the diagnostic. Used only for postprocessing purposes.
 every¶
Number of timesteps between each output or a time selection.
 flush_every¶
 Default
1
Number of timesteps or a time selection.
When
flush_every
coincides withevery
, the output file is actually written (“flushed” from the buffer). Flushing too often can dramatically slow down the simulation.
 time_average¶
 Default
1
(no averaging)
The number of timesteps for timeaveraging.
 fields¶
 Default
[]
(all fields are written)
List of the field names that are saved. By default, they all are. The full list of fields that are saved by this diagnostic:
BxByBzComponents of the magnetic fieldBx_mBy_mBz_mComponents of the magnetic field (timecentered)ExEyEzComponents of the electric fieldJxJyJzComponents of the total currentJx_abcJy_abcJz_abcComponents of the current due to species “abc”RhoRho_abcTotal densityDensity of species “abc”In
AMcylindrical
geometry, thex
,y
andz
indices are replaced byl
(longitudinal),r
(radial) andt
(theta). In addition, the angular Fourier modes are denoted by the suffix_mode_i
wherei
is the mode number. If a field is specified without its associated mode number, all available modes will be included. In summary, the list of fields reads as follows.Bl_mode_0, Bl_mode_1, etc.Br_mode_0, Br_mode_1, etc.Bt_mode_0, Bt_mode_1, etc.Components of the magnetic fieldEl_mode_0, El_mode_1, etc.Er_mode_0, Er_mode_1, etc.Et_mode_0, Et_mode_1, etc.Components of the electric fieldThe same notation works for Jl, Jr, Jt, and Rho
In the case of an envelope model for the laser (see Laser envelope model), the following fields are also available:
Env_A_absModule of laser vector potential’s complex envelope\(\tilde{A}\) (component along the transversedirection)Env_ChiTotal susceptibility \(\chi\)Env_E_absModule of laser electric field’s complex envelope\(\tilde{E}\) (component along the transversedirection)Env_Ex_absModule of laser electric field’s complex envelope\(\tilde{E}_x\) (component along the propagationdirection)
Note
To write these last three envelope fields with this diagnostics in "AMcylindrical"
geometry,
a dedicated block DiagFields
must be defined, e.g. with fields = ["Env_A_abs_mode_0", "Env_Chi_mode_0"]
.
 subgrid¶
 Default
None
(the whole grid is used)
A list of slices indicating a portion of the simulation grid to be written by this diagnostic. This list must have as many elements as the simulation dimension. For example, in a 3D simulation, the list has 3 elements. Each element can be:
None
, to select the whole grid along that dimensionan integer, to select only the corresponding cell index along that dimension
a python slice object to select regularlyspaced cell indices along that dimension.
This can be easily implemented using the numpy.s_ expression. For instance, in a 3D simulation, the following subgrid selects only every other element in each dimension:
from numpy import s_ DiagFields( #... subgrid = s_[::2, ::2, ::2] )
while this one selects cell indices included in a contiguous parallelepiped:
subgrid = s_[100:300, 300:500, 300:600]
Probe diagnostics¶
The fields from the previous section are taken at the PIC grid locations, but it is also possible to obtain the fields at arbitrary locations. These are called probes.
A probe interpolates the fields at either one point (0D), several points arranged in a line (1D), or several points arranged in a 2D or 3D grid.
Note
Probes follow the moving window. To obtain the fields at fixed points in the plasma instead, create a cold, chargeless species, and track the particles.
In “AMcylindrical” geometry, probes are defined with 3D Cartesian coordinates and cannot be separated per mode. Use Field diagnostics for cylindrical coordinates and information per mode.
To add one probe diagnostic, include the block DiagProbe
:
DiagProbe(
#name = "my_probe",
every = 10,
origin = [1., 1.],
corners = [
[1.,10.],
[10.,1.],
]
number = [100, 100],
fields = ["Ex", "Ey", "Ez"]
)
 name¶
Optional name of the diagnostic. Used only for postprocessing purposes.
 every¶
Number of timesteps between each output or a time selection.
 flush_every¶
 Default
1
Number of timesteps or a time selection.
When
flush_every
coincides withevery
, the output file is actually written (“flushed” from the buffer). Flushing too often can dramatically slow down the simulation.
 origin¶
 Type
A list of floats, of length equal to the simulation dimensionality.
The coordinates of the origin of the probe grid
 corners¶
 vectors¶
 Type
A list of lists of floats.
Defines the corners of the probe grid. Each corner is a list of coordinates (as many as the simulation dimensions).
When using
corners
, the absolute coordinates of each corner must be specified. When usingvectors
, the coordinates relative toorigin
must be specified.
 number¶
 Type
A list of integers, one for each dimension of the probe.
The number of points in each probe axis. Must not be defined for a 0D probe.
 fields¶
 Default
[]
, which means["Ex", "Ey", "Ez", "Bx", "By", "Bz", "Jx", "Jy", "Jz", "Rho"]
A list of fields among:
the electric field components
"Ex"
,"Ey"
,"Ez"
the magnetic field components
"Bx"
,"By"
,"Bz"
the Poynting vector components
"PoyX"
,"PoyY"
,"PoyZ"
the current density components
"Jx"
,"Jy"
,"Jz"
and density"Rho"
the current density
"Jx_abc"
,"Jy_abc"
,"Jz_abc"
and density"Rho_abc"
of a given species named"abc"
In the case of an envelope model for the laser (see Laser envelope model), the following fields are also available:
"Env_A_abs"
,"Env_Chi"
,"Env_E_abs"
,"Env_Ex_abs"
.
 time_integral¶
 Default
False
If
True
, the output is integrated over time. As this option forces field interpolation at every timestep, it is recommended to use few probe points.
Examples of probe diagnostics
0D probe in 1D simulation
DiagProbe( every = 1, origin = [1.2] )
1D probe in 1D simulation
DiagProbe( every = 1, origin = [1.2], corners = [[5.6]], number = [100] )
1D probe in 2D simulation
DiagProbe( every = 1, origin = [1.2, 4.], corners = [[5.6, 4.]], number = [100] )
2D probe in 2D simulation
DiagProbe( every = 1, origin = [0., 0.], corners = [ [10.,0.], [0.,10.] ], number = [100, 100] )
ParticleBinning diagnostics¶
A particle binning diagnostic collects data from the macroparticles and processes them during runtime. It does not provide information on individual particles: instead, it produces averaged quantities like the particle density, currents, etc.
The data is discretized inside a “grid” chosen by the user. This grid may be of any dimension.
Examples:
1dimensional grid along the position \(x\) (gives density variation along \(x\))
2dimensional grid along positions \(x\) and \(y\) (gives density map)
1dimensional grid along the velocity \(v_x\) (gives the velocity distribution)
2dimensional grid along position \(x\) and momentum \(p_x\) (gives the phasespace)
1dimensional grid along the kinetic energy \(E_\mathrm{kin}\) (gives the energy distribution)
3dimensional grid along \(x\), \(y\) and \(E_\mathrm{kin}\) (gives the density map for several energies)
1dimensional grid along the charge \(Z^\star\) (gives the charge distribution)
0dimensional grid (simply gives the total integrated particle density)
Each dimension of the grid is called “axis”.
You can add a particle binning diagnostic by including a block DiagParticleBinning()
in the namelist,
for instance:
DiagParticleBinning(
#name = "my binning",
deposited_quantity = "weight",
every = 5,
time_average = 1,
species = ["electrons1", "electrons2"],
axes = [
["x", 0., 10, 100],
["ekin", 0.1, 100, 1000, "logscale", "edge_inclusive"]
]
)
 name¶
Optional name of the diagnostic. Used only for postprocessing purposes.
 deposited_quantity¶
The type of data that is summed in each cell of the grid. Consider reading this to understand the meaning of the
weight
."weight"
results in a number density."weight_charge"
results in a charge density."weight_charge_vx"
results in the \(j_x\) current density (same with \(y\) and \(z\))."weight_p"
results in the momentum density (same with \(p_x\), \(p_y\) and \(p_z\))."weight_ekin"
results in the energy density."weight_vx_px"
results in thexx
pressure (same with yy, zz, xy, yz and xz)."weight_chi"
results in the quantum parameter density (only for species with radiation losses).with a userdefined python function, an arbitrary quantity can be calculated (the numpy module is necessary). This function should take one argument, for instance
particles
, which contains the attributesx
,y
,z
,px
,py
,pz
,charge
,weight
,chi
andid
. Each of these attributes is a numpy array containing the data of all particles in one patch. The function must return a numpy array of the same shape, containing the desired deposition of each particle. For example, defining the following function:def stuff(particles): return particles.weight * particles.px
passed as
deposited_quantity=stuff
, the diagnostic will sum the weights \(\times\; p_x\).You may also pass directly an implicit (lambda) function using:
deposited_quantity = lambda p: p.weight * p.px
 every¶
The number of timesteps between each output, or a time selection.
 flush_every¶
 Default
1
Number of timesteps or a time selection.
When
flush_every
coincides withevery
, the output file is actually written (“flushed” from the buffer). Flushing too often can dramatically slow down the simulation.
 time_average¶
 Default
1
The number of timesteps during which the data is averaged before output.
 species¶
A list of one or several species’
name
. All these species are combined into the same diagnostic.
 axes¶
A list of axes that define the grid. There may be as many axes as wanted (there may be zero axes).
Syntax of one axis:
[type, min, max, nsteps, "logscale", "edge_inclusive"]
type
is one of:"x"
,"y"
,"z"
: spatial coordinates ("moving_x"
with a moving window)"px"
,"py"
,"pz"
,"p"
: momenta"vx"
,"vy"
,"vz"
,"v"
: velocities"gamma"
,"ekin"
: energies"chi"
: quantum parameter"charge"
: the particles’ electric chargeor a python function with the same syntax as the
deposited_quantity
. Namely, this function must accept one argument only, for instanceparticles
, which holds the attributesx
,y
,z
,px
,py
,pz
,charge
,weight
andid
. Each of these attributes is a numpy array containing the data of all particles in one patch. The function must return a numpy array of the same shape, containing the desired quantity of each particle that will decide its location in the histogram binning.
The axis is discretized for
type
frommin
tomax
innsteps
bins.The
min
andmax
may be set to"auto"
so that they are automatically computed from all the particles in the simulation. This option can be bad for performances.The optional keyword
logscale
sets the axis scale to logarithmic instead of linear (bins become uneven).The optional keyword
edge_inclusive
includes the particles outside the range [min
,max
] into the extrema bins.
Examples of particle binning diagnostics
Variation of the density of species
electron1
from \(x=0\) to 1, every 5 timesteps, without timeaveragingDiagParticleBinning( deposited_quantity = "weight", every = 5, time_average = 1, species = ["electron1"], axes = [ ["x", 0., 1., 30] ] )
Density map from \(x=0\) to 1, \(y=0\) to 1
DiagParticleBinning( deposited_quantity = "weight", every = 5, time_average = 1, species = ["electron1"], axes = [ ["x", 0., 1., 30], ["y", 0., 1., 30] ] )
Velocity distribution from \(v_x = 0.1\) to \(0.1\)
DiagParticleBinning( deposited_quantity = "weight", every = 5, time_average = 1, species = ["electron1"], axes = [ ["vx", 0.1, 0.1, 100] ] )
Phase space from \(x=0\) to 1 and from \(px=1\) to 1
DiagParticleBinning( deposited_quantity = "weight", every = 5, time_average = 1, species = ["electron1"], axes = [ ["x", 0., 1., 30], ["px", 1., 1., 100] ] )
Energy distribution from 0.01 to 1 MeV in logarithmic scale. Note that the input units are \(m_ec^2 \sim 0.5\) MeV
DiagParticleBinning( deposited_quantity = "weight", every = 5, time_average = 1, species = ["electron1"], axes = [ ["ekin", 0.02, 2., 100, "logscale"] ] )
\(x\)\(y\) density maps for three bands of energy: \([0,1]\), \([1,2]\), \([2,\infty]\). Note the use of
edge_inclusive
to reach energies up to \(\infty\)DiagParticleBinning( deposited_quantity = "weight", every = 5, time_average = 1, species = ["electron1"], axes = [ ["x", 0., 1., 30], ["y", 0., 1., 30], ["ekin", 0., 6., 3, "edge_inclusive"] ] )
Charge distribution from \(Z^\star =0\) to 10
DiagParticleBinning( deposited_quantity = "weight", every = 5, time_average = 1, species = ["electron1"], axes = [ ["charge", 0.5, 10.5, 11] ] )
Screen diagnostics¶
A screen collects data from the macroparticles when they cross a surface. It processes this data similarly to the particle binning diagnostics as it makes a histogram of the macroparticle properties. The only difference is that the histogram is made only by the particles that cross the surface.
You can add a screen by including a block DiagScreen()
in the namelist,
for instance:
DiagScreen(
#name = "my screen",
shape = "plane",
point = [5., 10.],
vector = [1., 0.],
direction = "canceling",
deposited_quantity = "weight",
species = ["electron"],
axes = [["a", 10.*l0, 10.*l0, 40],
["px", 0., 3., 30]],
every = 10
)
 name¶
Optional name of the diagnostic. Used only for postprocessing purposes.
 shape¶
The shape of the screen surface:
"plane"
,"sphere"
, or"cylinder"
.
 point¶
 Type
A list of floats
[X]
in 1D,[X,Y]
in 2D,[X,Y,Z]
in 3D
The coordinates of a point that defines the screen surface: a point of the
"plane"
, the center of the"sphere"
, or a point on the"cylinder"
axis.
 vector¶
 Type
A list of floats
[X]
in 1D,[X,Y]
in 2D,[X,Y,Z]
in 3D
The coordinates of a vector that defines the screen surface: the normal to the
"plane"
, a radius of the"sphere"
. or the axis of the"cylinder"
(in the latter case, the vector norm defines the cylinder radius).
 direction¶
 Default
"both"
Determines how particles are counted depending on which side of the screen they come from.
"both"
to account for both sides."forward"
for only the ones in the direction of thevector
."backward"
for only the ones in the opposite direction."canceling"
to count negatively the ones in the opposite direction.
 deposited_quantity¶
Identical to the
deposited_quantity
of particle binning diagnostics.
 every¶
The number of timesteps between each output, or a time selection.
 flush_every¶
 Default
1
Number of timesteps or a time selection.
When
flush_every
coincides withevery
, the output file is actually written (“flushed” from the buffer). Flushing too often can dramatically slow down the simulation.
 species¶
A list of one or several species’
name
. All these species are combined into the same diagnostic.
 axes¶
A list of “axes” that define the grid of the histogram. It is identical to that of particle binning diagnostics, with the addition of four types of axes:
If
shape="plane"
, then"a"
and"b"
are the axes perpendicular to thevector
.If
shape="sphere"
, then"theta"
and"phi"
are the angles with respect to thevector
.If
shape="cylinder"
, then"a"
is along the cylinder axis and"phi"
is the angle around it.
RadiationSpectrum diagnostics¶
A radiation spectrum diagnostic computes (at a given time) the instantaneous power spectrum following from the incoherent emission of highenergy photons by accelerated charge (see Highenergy photon emission & radiation reaction for more details on the emission process and its implementation in Smilei).
It is similar to the particle binning diagnostics, with an extra axis of binning: the emitted photon energy. The other axes remain available to the user.
A radiation spectrum diagnostic is defined by a block RadiationSpectrum()
:
DiagRadiationSpectrum(
#name = "my radiation spectrum",
every = 5,
flush_every = 1,
time_average = 1,
species = ["electrons1", "electrons2"],
photon_energy_axis = [0., 1000., 100, 'logscale'],
axes = []
)
 name¶
Optional name of the diagnostic. Used only for postprocessing purposes.
 every¶
The number of timesteps between each output, or a time selection.
 flush_every¶
 Default
1
Number of timesteps or a time selection.
When
flush_every
coincides withevery
, the output file is actually written (“flushed” from the buffer). Flushing too often can dramatically slow down the simulation.
 time_average¶
 Default
1
The number of timesteps during which the data is averaged before output.
 species¶
A list of one or several species’
name
that emit the radiation. All these species are combined into the same diagnostic.
 photon_energy_axis¶
The axis of photon energies (in units of \(m_e c^2\)). The syntax is similar to that of particle binning diagnostics.
Syntax:
[min, max, nsteps, "logscale"]
 axes¶
An additional list of “axes” that define the grid. There may be as many axes as wanted (there may be zero axes). Their syntax is the same that for “axes” of a particle binning diagnostics.
Examples of radiation spectrum diagnostics
Timeintegrated over the full duration of the simulation:
DiagRadiationSpectrum( every = Nt, time_average = Nt, species = ["electrons"], photon_energy_axis = [0., 1000., 100, 'logscale'], axes = [] )
Angularlyresolved instantaneous radiation spectrum. The diagnostic considers that all electrons emit radiation in the direction of their velocity:
from numpy import arctan2, pi def angle(p): return arctan2(p.py,p.px) DiagRadiationSpectrum( every = 10, species = ["electrons"], photon_energy_axis = [0., 1000., 100, 'logscale'], axes = [ [angle,pi,pi,90] ] )
TrackParticles diagnostics¶
A particle tracking diagnostic records the macroparticle positions and momenta at various timesteps. Typically, this is used for plotting trajectories.
You can add a tracking diagnostic by including a block DiagTrackParticles()
in the namelist,
for instance:
DiagTrackParticles(
species = "electron",
every = 10,
# flush_every = 100,
# filter = my_filter,
# attributes = ["x", "px", "py", "Ex", "Ey", "Bz"]
)
 every¶
 Default
0
Number of timesteps between each output of particles trajectories, or a time selection. If nonzero, the particles positions will be tracked and written in a file named
TrackParticlesDisordered_abc.h5
(whereabc
is the species’name
).
 flush_every¶
 Default
1
Number of timesteps or a time selection.
When
flush_every
coincides withevery
, the output file for tracked particles is actually written (“flushed” from the buffer). Flushing too often can dramatically slow down the simulation.
 filter¶
A python function giving some condition on which particles are tracked. If none provided, all particles are tracked. To use this option, the numpy package must be available in your python installation.
The function must have one argument, that you may call, for instance,
particles
. This object has several attributesx
,y
,z
,px
,py
,pz
,charge
,weight
andid
. Each of these attributes are provided as numpy arrays where each cell corresponds to one particle. The function must return a boolean numpy array of the same shape, containingTrue
for particles that should be tracked, andFalse
otherwise.The following example selects all the particles that verify \(1<p_x<1\) or \(p_z>3\):
def my_filter(particles): return (particles.px>1.)*(particles.px<1.) + (particles.pz>3.)
Warning
The px
, py
and pz
quantities are not exactly the momenta.
They are actually the velocities multiplied by the lorentz factor, i.e.,
\(\gamma v_x\), \(\gamma v_y\) and \(\gamma v_z\). This is true only
inside the filter
function (not for the output of the diagnostic).
Note
The id
attribute contains the particles identification number.
This number is set to 0 at the beginning of the simulation. Only after particles have
passed the filter, they acquire a positive id
.
Note
For advanced filtration, Smilei provides the quantity Main.iteration
,
accessible within the filter
function. Its value is always equal to the current
iteration number of the PIC loop. The current time of the simulation is thus
Main.iteration * Main.timestep
.
 attributes¶
 Default
["x","y","z","px","py","pz","w"]
A list of strings indicating the particle attributes to be written in the output. The attributes may be the particles’ spatial coordinates (
"x"
,"y"
,"z"
), their momenta ("px"
,"py"
,"pz"
), their electrical charge ("q"
), their statistical weight ("w"
), their quantum parameter ("chi"
, only for species with radiation losses) or the fields interpolated at their positions ("Ex"
,"Ey"
,"Ez"
,"Bx"
,"By"
,"Bz"
).
Performances diagnostics¶
The performances diagnostic records information on the computational load and timers for each MPI process or for each patch in the simulation.
Only one block DiagPerformances()
may be added in the namelist, for instance:
DiagPerformances(
every = 100,
# flush_every = 100,
# patch_information = True,
)
 every¶
 Default
0
Number of timesteps between each output, or a time selection.
 flush_every¶
 Default
1
Number of timesteps or a time selection.
When
flush_every
coincides withevery
, the output file is actually written (“flushed” from the buffer). Flushing too often might dramatically slow down the simulation.
 patch_information¶
 Default
False
If
True
, some information is calculated at the patch level (seePerformances()
) but this may impact the code performances.
Time selections¶
Several components (mainly diagnostics) may require a selection of timesteps to
be chosen by the user. When one of these timesteps is reached, the diagnostics will
output data. A time selection is given through the parameter every
and is a list
of several numbers.
You may chose between five different syntaxes:
every = [ period ] # Syntax 1
every = [ start, period ] # Syntax 2
every = [ start, end, period ] # Syntax 3
every = [ start, end, period, repeat ] # Syntax 4
every = [ start, end, period, repeat, spacing ] # Syntax 5
where
start
is the first timestep of the selection (defaults to 0);end
is the last timestep of the selection (defaults to ∞);period
is the separation between outputs (defaults to 1);repeat
indicates how many outputs to do at each period (defaults to 1);spacing
is the separation between each repeat (defaults to 1).
For more clarity, this graph illustrates the five syntaxes for time selections:
Tips
The syntax
every = period
is also accepted.Any value set to
0
will be replaced by the default value.Special case:
every=0
means no output.The numbers may be nonintegers (apart from
repeat
). The closest timesteps are chosen.
Checkpoints¶
The simulation state can be saved (dumped) at given times (checkpoints) in order to be later restarted at that point.
A few things are important to know when you need dumps and restarts.
Do not restart the simulation in the same directory as the previous one. Files will be overwritten, and errors may occur. Create a new directory for your restarted simulation.
Manage your disk space: each MPI process dumps one file, and the total can be significant.
The restarted runs must have the same namelist as the initial simulation, except the Checkpoints block, which can be modified.
Checkpoints(
# restart_dir = "dump1",
dump_step = 10000,
dump_minutes = 240.,
exit_after_dump = True,
keep_n_dumps = 2,
)
Parameters to save the state of the current simulation
 dump_step¶
 Default
0
The number of timesteps between each dump. If
0
, no dump is done.
 dump_minutes¶
 Default
0.
The number of minutes between each dump. If
0.
, no dump is done.May be used in combination with
dump_step
.
 exit_after_dump¶
 Default
True
If
True
, the code stops after the first dump. IfFalse
, the simulation continues.
 keep_n_dumps¶
 Default
2
This tells Smilei to keep, in the current run, only the last
n
dumps. Older dumps will be overwritten.The default value,
2
, saves one extra dump in case of a crash during the next dump.
 file_grouping¶
 Default
0
(no grouping)The maximum number of checkpoint files that can be stored in one directory. Subdirectories are created to accomodate for all files. This is useful on filesystem with a limited number of files per directory.
 dump_deflate¶
to do
Parameters to restart from a previous simulation
 restart_dir¶
 Default
None
The directory of a previous run from which Smilei should restart. For the first run, do not specify this parameter.
This path must either absolute or be relative to the current directory.
Note
In many situations, the restarted runs will have the exact same namelist as the initial simulation, except this
restart_dir
parameter, which points to the previous simulation folder. You can use the same namelist file, and simply add an extra argument when you launch the restart:
mpirun ... ./smilei mynamelist.py "Checkpoints.restart_dir='/path/to/previous/run'"
 restart_number¶
 Default
None
The number of the dump (in the previous run) that should be used for the restart. For the first run, do not specify this parameter.
In a previous run, the simulation state may have been dumped several times. These dumps are numbered 0, 1, 2, etc. until the number
keep_n_dumps
.
Variables defined by Smilei¶
Smilei passes the following variables to the python interpreter for use in the namelist. They should not be redefined by the user!
 smilei_mpi_rank¶
The MPI rank of the current process.
 smilei_mpi_size¶
The total number of MPI processes.