PaperlayerView exact revision on arXiv

Source-linked AI summary

SMILEI: a collaborative, open-source, multi-purpose particle-in-cell code for plasma simulation

J. Derouillat, A. Beck, F. Pérez, T. Vinci, M. Chiaramello, A. Grassi, M. Flé, G. Bouchard, I. Plotnikov, N. Aunai, J. Dargent, C. Riconda, M. Grech

arXiv:1702.05128v1physics.plasm-ph

TL;DR

Plasma simulations must support diverse physics while exploiting rapidly evolving massively parallel HPC architectures. Smilei addresses this need with a collaborative, open-source PIC code combining multi-purpose physics modules with patch-based hybrid parallelization. The code scales to several hundred thousand computing elements and has been tested across multiple supercomputers, architectures, and processors.

  • Problem

    Plasma codes face difficulties adapting to massively parallel and emerging exascale architectures while supporting applications from laser-plasma interaction to astrophysics.

  • Method

    Smilei combines a collaborative physicist-HPC development model with open-source PIC physics, patch-based MPI-OpenMP parallelization, compact patch distribution, and dynamic load balancing.

  • Results

    Smilei scales well to several hundred thousand computing elements and was tested on multiple supercomputers, architectures, and processors.

  • Takeaways & Limitations

    Smilei provides a versatile, multi-purpose plasma-simulation tool with an object-oriented C++ structure and user-friendly Python interface.

  • Takeaways & Limitations

    Simulations using ionization, binary collisions, or collisional ionization cannot be arbitrarily scaled because their reference time normalization must be defined in SI units.

Abstract

from arXiv · show

SMILEI is a collaborative, open-source, object-oriented (C++) particle-in-cell code. To benefit from the latest advances in high-performance computing (HPC), SMILEI is co-developed by both physicists and HPC experts. The code's structures, capabilities, parallelization strategy and performances are discussed. Additional modules (e.g. to treat ionization or collisions), benchmarks and physics highlights are also presented. Multi-purpose and evolutive, SMILEI is applied today to a wide range of physics studies, from relativistic laser-plasma interaction to astrophysical plasmas.

Program summary

Smilei simulates plasma kinetics with the PIC method and supports additional physics modules. Its hybrid MPI-OpenMP, patch-based parallelization targets efficient execution on massively parallel supercomputers.

  • Plasma kinetic simulation supports studies ranging from laser-plasma interaction to astrophysics.
  • The Vlasov-Maxwell system for collisionless plasmas is solved using the Particle-In-Cell method.Additional modules account for effects such as collisions and ionization.
  • A hybrid MPI-OpenMP strategy based on patch-based super-decomposition improves cache use, dynamic load balancing, and performance on massively parallel supercomputers.

1. Introduction

Smilei was created as an open-source, collaborative PIC code for increasingly demanding laser-plasma and broader plasma simulations. Its development combines physicists’ domain knowledge with HPC expertise to address evolving supercomputer architectures.

  • PIC methods are central to kinetic plasma simulations across laboratory, accelerator, semiconductor, cosmological, and astrophysical studies.
  • Smilei is an open-source PIC code developed collaboratively by physicists and HPC experts for new HPC architectures.
  • Ultra-high-intensity laser facilities motivate Smilei because laser-plasma interaction enables particle acceleration, ultrashort bright light sources, and pair production.
  • The rapid growth of massively parallel supercomputer cores challenges existing plasma codes, requiring stronger collaboration between physicists and HPC specialists.
  • Smilei is intended as a multi-purpose, collaborative code spanning laser-plasma interaction to astrophysics.
  • The paper surveys Smilei’s principles, structure, performance, capabilities, benchmarks, examples, parallelization, additional modules, and input interface.

2. The Particle-In-Cell (PIC) method for collisionless plasmas

The PIC method represents plasma distributions with weighted quasi-particles, advances them under electromagnetic forces, and solves Maxwell’s equations on a discrete grid. Smilei uses normalized units and supports collisionless-plasma modeling with explicitly added collision modules.

  • 2.1. The Maxwell-Vlasov model: The collisionless-plasma model uses the Vlasov-Maxwell system, where particle distributions evolve under the Lorentz force and modify collective electromagnetic fields.
  • 2.2. Reference units: Smilei’s fully relativistic PIC formulation normalizes velocity to c, charge to e, mass to m_e, momentum to m_ec, and energy to m_ec^2.
  • 2.2. Reference units: Time and space use a reference frequency ω_r, with distances normalized to c/ω_r and fields and densities assigned corresponding reference units.
  • 2.1. The Maxwell-Vlasov model: Standard PIC uses macroscopic electromagnetic fields and therefore does not include particle collisions without an additional ad hoc module.
  • 2.3. Quasi-particles and the PIC method: The PIC method discretizes each distribution function as weighted quasi-particles with positions, momenta, Dirac distributions, and shape functions.
  • 2.3. Quasi-particles and the PIC method: Particle trajectories are integrated continuously while Maxwell’s equations are solved on a discrete spatial grid divided into cells.
  • 2.3. Quasi-particles and the PIC method: Ghost cells are required at domain boundaries so finite-size quasi-particles’ charge and current densities project correctly onto the grid.

2.4. Time- and space-centered discretization

Smilei discretizes electromagnetic fields on a staggered Yee grid and advances particles with leap-frog time centering. This arrangement centers spatial curl operators and aligns field values with particle-force calculations.

  • The FDTD approach and refined variants discretize electromagnetic fields on a staggered Yee grid.
  • The Yee grid spatially centers the discretized curl operators in Maxwell’s equations.
  • Table 2 summarizes Smilei’s PIC algorithm.
  • Charge and current densities are located according to the convention of defining charge densities at cell nodes.
  • Leap-frog time advancement places particle positions at integer time steps and velocities at half-integer time steps.

2.5. Initialization of the simulation

Smilei initializes PIC simulations by loading weighted particles, projecting initial charge and current densities, and computing electromagnetic fields across the grid. Users can define these conditions directly or restart from a previous simulation.

  • Initialization workflow: Initialization loads particles, computes initial charge and current densities on the grid, and computes initial electric and magnetic fields.All three steps can be performed from a user-defined input file or by restarting from a previous simulation.
  • Particle loading: Users prescribe each species’ density, particles per cell, mean velocity, and temperature before particle loading.Positions are uniformly distributed within each cell, while momenta are sampled from a requested distribution.
  • Particle loading: Smilei supports zero-temperature, Maxwellian, and Maxwell-Jüttner momentum distributions, including relativistic drift loading for the latter.The Maxwell-Jüttner distribution is the relativistic generalization of the Maxwellian distribution.
  • Particle weighting: Each particle receives a numerical weight based on the density of its originating cell, benefiting initially highly inhomogeneous density distributions.Variable weighting adapts particle representation to the cell-associated density.
  • Field initialization: Initial charge and current densities are projected onto grid points from the charge and current contained in surrounding cells.The direct projection uses the cell geometry and spatial dimensions of the simulation.
  • Field initialization: Initial electric fields are obtained by solving Poisson’s equation with a conjugate-gradient method suited to massively parallel computers.The method mainly requires local information exchange between adjacent computational domains.
  • Field initialization: External divergence-free electric or magnetic fields may be added if they satisfy Maxwell’s equations, including Gauss’ and Poisson’s equations.

2.6. The PIC loop

The PIC loop advances particles and fields through repeated interpolation, particle pushing, charge-conserving projection, and Maxwell updates. Smilei supports alternative particle pushers, field solvers, boundary conditions, and discretization choices within this cycle.

  • Loop sequence: Each time-step interpolates fields to particles, updates particle velocities and positions, projects charge and currents, and advances electromagnetic fields.
  • Particle pusher: Smilei computes particle momentum and position updates with a second-order leap-frog integrator using Boris or Vay schemes.Both schemes are implemented for quasi-particles after electromagnetic fields are known at their positions.
  • Particle pusher: Vay’s solver addresses relativistic orbit errors that Boris can introduce in special electromagnetic configurations, with an additional but not large computational cost.
  • Charge deposition: Esirkepov’s charge-conserving algorithm computes grid-direction current densities from particle charge flux through cell borders.The particle current in the z-direction is computed by direct projection when z is not a resolved grid dimension.
  • Charge deposition: Charge density may be projected for diagnostics, while charge-conserving field advancement does not use it to advance electromagnetic fields.
  • Field solver: Maxwell-Ampère and Maxwell-Faraday updates advance electric and magnetic fields after the currents are known.The fields are discretized on a staggered Yee grid, with extended curl stencils available for Faraday derivatives.
  • Field solver: Charge-conserving deposition preserves the divergence constraints when the initial Gauss and Poisson equations are satisfied.This property motivates using Esirkepov’s projection scheme.
  • Field solver: FDTD solvers are subject to a Courant-Friedrich-Lewy condition requiring the time-step to remain below a solver-dependent limit.The standard FDTD solver corresponds to αx = 1 and ηx = βxy = βxz = 0.

3. An evolutive, multi-purpose code

Smilei uses object-oriented C++ design and polymorphism to separate physics, data structures, computing algorithms, and parallelism. This supports one evolutive executable across dimensions and physics while addressing HPC data-management demands.

  • Design objectives: Smilei’s C++ approach separates physics from computing aspects to support performance, maintainability, a large user community, and varied applications.
  • C++ architecture: Object-oriented C++ structures components independently and avoids large parameter lists and global variables that are inefficient in parallel computing.
  • C++ architecture: Particle species and electromagnetics classes organize quasi-particle data, fields, and operations such as pushers and boundary conditions.Interpolators transfer field data to particle positions, while projectors transfer particle data to grid points.
  • Polymorphism: Polymorphic field classes such as Field1D and Field2D share inherited functions while using dimension-specific data structures selected through a factory design pattern.
  • Polymorphism: A single executable supports different dimensions, interpolation orders, and physics components without requiring many code versions.
  • Polymorphism: Excessive virtualization or many object layers can impose significant computational cost, so Smilei passes whole data structures to operators before casting them as needed.
  • HPC adaptability: Wrappers decouple computing algorithms from raw data structures, while separate parallelism levels decouple parallelism management from physics calculations.
  • I/O management: Parallel I/O addresses communication, memory, and post-processing overheads that grow with simulation size; Smilei’s preliminary benchmark reached 2.6 Gb/s.The benchmark wrote 60 Gb in 24 seconds per selected timestep for a 30720 x 15360-cell domain.

4. Parallelization

SMILEI targets exascale-era plasma simulation by combining fine-grained patch decomposition with hybrid MPI–OpenMP parallelism, dynamic load balancing, and cache-aware data locality. The strategy scales well for homogeneous plasmas, while hybrid balancing mitigates—but does not eliminate—load imbalance in evolving laser-plasma simulations.

  • Parallelization requirements: Exascale-oriented parallelization must expose massive parallelism while preserving data locality and minimizing load imbalance.SMILEI addresses this through a patch-based strategy developed for evolving HPC architectures.
  • Patch-based decomposition: Particle-based decomposition offers load balance but causes particle communications and random grid access, whereas grid decomposition complicates balance as particles move.SMILEI therefore uses fine-grain patch decomposition, storing particles with the grid portions they interact with.
  • Hybrid MPI–OpenMP: Hybrid MPI–OpenMP patching provides scalability and dynamic load balancing, while reducing global communication through fewer, larger MPI domains.OpenMP dynamically schedules patches within each MPI process; inter-process synchronization remains more expensive than shared-memory synchronization.
  • Pure MPI: MPI performs well for homogeneous plasmas when MPI domains are not too small relative to the L1 cache.This condition defines the practical boundary of the pure-MPI approach in the reported strong-scaling study.
  • Hybrid performance: A nearly 20-fold wall-clock increase under pure MPI occurs after load imbalance develops, while the best 64×12 hybrid case still loses more than a factor of 4.OpenMP reduces sensitivity to imbalance, but cannot balance work across MPI domains; increasing threads also incurs synchronization costs across processors.
  • Dynamic load balancing: Dynamic load balancing adds almost 40% gain over the previous best 64×12 case and delivers similar results with only 6 OpenMP threads.More patches improve balancing under strong imbalance but increase synchronization and ghost-cell costs when imbalance is weak.

5. Additional modules

Smilei extends the standard collisionless PIC model with user-selectable modules for numerical stability, ionization, collisions, and external currents. These additions support relativistic and strongly coupled plasma simulations while imposing specific modeling and scaling constraints.

  • Overview: Smilei provides additional modules to address user needs, including numerical filters, antennas, field ionization, and collisional processes.The modules extend the core PIC treatment beyond collisionless plasma dynamics.
  • Numerical stability: FDTD-based PIC methods can produce grid-Cerenkov instability for relativistic beams or drifting plasmas, so Smilei implements electric-field and current-density filtering.The instability arises from numerical dispersion that reduces the computed light-wave velocity.
  • Numerical stability: Both filtering methods, used separately or together, satisfactorily reduce grid-Cerenkov instability in relativistically drifting electron-positron plasmas.The methods were applied in collisionless shock studies.
  • Ionization: Smilei treats field ionization with a Monte-Carlo module based on quasi-static tunnel-ionization rates and supports multiple ionization events within one timestep.The implementation computes the rate using the |m| = 0 channel and applies ionization after field interpolation and before particle pushing.
  • Limitations: Simulations involving field ionization, binary collisions, or collisional ionization cannot be arbitrarily scaled because their reference time normalization is constrained.The paper explicitly identifies this as a limitation of simulations using these additional physics modules.
  • Collisions and collisional ionization: Smilei includes relativistic collision modeling and reports stopping-power agreement with theoretical calculations for high-energy electrons traversing ionized hydrogen plasma.The collision model incorporates relativistic particles, low-temperature corrections, and a variable Coulomb logarithm.
  • Collisions and collisional ionization: The multiple-ionization scheme matches well-resolved simulations better than single-ionization at disparate timesteps, achieving comparable accuracy without an order-of-magnitude timestep reduction.The scheme was tested across materials and incident electron energies, including zinc-plasma simulations.

6. User interface

Smilei combines Python-based namelists with runtime diagnostics to make simulation setup, data collection, and analysis configurable within the code. Users can initialize complex simulations and selectively record global, field, probe, trajectory, and distribution data.

  • Python interface: Smilei’s Python namelist lets users execute complex initialization operations and generate arbitrary simulation elements at runtime.The interpreter remains active during initialization and can remain active when functions are evaluated throughout the simulation.
  • Diagnostics: Diagnostics are integrated into the main code, execute at runtime, and can be scheduled at user-defined simulation times.They are not separate post-processing modules.
  • Diagnostics: Scalar diagnostics combine field and particle data across processors and write quantities such as energy balance, averaged particle properties, and global field extrema.The output is stored in a dedicated scalar file.
  • Diagnostics: Field diagnostics output concatenated electromagnetic, density, and current arrays, with optional temporal averaging over requested timesteps.These diagnostics provide direct copies of code arrays after processor concatenation.
  • Diagnostics: Probe diagnostics reduce storage by recording interpolated fields at user-defined points arranged as lines, rectangles, or parallelepipeds.They avoid writing whole simulation arrays when localized sampling is sufficient.
  • Diagnostics: Tracking stores independently configurable histories for selected particle species, assigning each tracked particle a unique persistent number.Custom output frequencies can be specified per species.
  • Diagnostics: Particle diagnostics deposit quasi-particle weights onto user-defined multidimensional grids to produce distributions and quantities such as current or pressure components.These diagnostics also support temporal averaging.

7. Physics highlights

Smilei is demonstrated across laser-plasma and astrophysical simulations, covering harmonic generation, pulse amplification, magnetopause reconnection, and collisionless shocks. The examples report both application-specific observables and plasma-kinetic behavior.

  • Laser-plasma interaction: A two-dimensional laser-solid simulation produces high harmonics up to order 16 and ejects electrons with energies up to 10 MeV.The energetic electrons predominantly propagate along directions close to the reflected laser direction.
  • Laser-plasma interaction: In plasma amplification simulations, the seed reaches nearly 5× its initial intensity while largely conserving spot size and phase front.The reported final seed intensity is Iout_s ≃4.6×10^15 W/cm2.
  • Laser-plasma interaction: A second amplification configuration reaches Iout_s ≃3×10^15 W/cm2, three times the initial intensity, while preserving a spot size of about 28 µm.The setup uses two pump lasers angled at ±6° to the x-axis.
  • Magnetic reconnection: A plasmaspheric plume reaches the magnetopause reconnection site at about t = 300 Ω_ci^-1 and subsequently fills the exhaust with mixed plume and solar-wind plasma.The transition includes important fluctuations caused by plasmoid formation.
  • Magnetic reconnection: When the plume reaches the reconnection site, the reconnection rate collapses to about half its previous value before a quasi-steady plume-reconnection regime develops.The transition phase lasts about 100 Ω_ci^-1.
  • Collisionless shocks: In a Weibel-mediated pair-plasma shock, magnetic turbulence stops incoming particles and produces downstream compression up to nd ≃3.2 n0, consistent with Rankine-Hugoniot predictions.The simulation develops filamentary structures ahead of the shock and a turbulent magnetic field at its front.
  • Collisionless shocks: The shocked downstream distribution initially becomes isotropic and quasithermal near Td ≃4.5 mec2, then develops a supra-thermal γ^-2.5 tail.The tail is identified as characteristic of first-order Fermi acceleration at the shock front.

8. Conclusions

Smilei is presented as a versatile, open-source PIC code combining an object-oriented C++ core, Python interface, modular physics, and scalable patch-based parallelization. Its applications span laser-plasma interaction and astrophysics, while development continues toward broader physics and software capabilities.

  • Contributions: Smilei’s object-oriented C++ structure and Python interface support a versatile, multi-purpose plasma simulation tool.The conclusion emphasizes the combination of maintainable code structure and user-friendly input.
  • Performance: Patch-based super-decomposition improves cache use, enables straightforward dynamic load balancing, and scales to several hundred thousand computing elements.The code was tested across Bullx and BlueGene/Q systems and multiple Intel and IBM processors.
  • Capabilities and use: Smilei includes Monte-Carlo modules for binary collisions and collisional and field ionization, and is used by a growing community.The presented applications cover both laser-plasma interaction and astrophysics.
  • Future development: Ongoing development targets quantum-electrodynamics modules, spectral Maxwell solvers, vectorization, documentation, and automated continuous-integration benchmarks.These efforts extend both the physics coverage and maintainability of the code.

Appendix A. Quasi-particles shape functions and interpolation/projection order

SMILEI represents quasi-particles with symmetric, normalized shape functions whose support spans a cell-based volume. The interpolation/projection order n determines the shape-function construction, from point-like particles to recursively defined higher-order functions.

  • The shape function is symmetric, normalized, and non-zero around x = 0 over a region extending n ∆xµ in each spatial direction.Its support defines the quasi-particle volume Vp = Πµn∆xµ.
  • The integer n is the interpolation/projection order used to characterize the quasi-particle shape.
  • Shape functions are constructed as products of one-dimensional functions over the D spatial dimensions.
  • For n = 0, the one-dimensional shape function is point-like and equals the Dirac distribution δ(x).
  • Higher-order shape functions with n > 0 are generated recursively, with explicit forms given up to order 4.

Field interpolation at the particle position

Field interpolation reconstructs the electric or magnetic field from grid values and evaluates it at each quasi-particle’s position. The same construction extends directly from one-dimensional grids to arbitrary spatial dimensions.

  • A quasi-particle at position xp sees an electric or magnetic field interpolated from the surrounding grid.
  • The grid field F(x) is reconstructed from values associated with grid points xi.
  • Substituting the grid reconstruction and the recursive shape function yields a direct interpolation formula at the quasi-particle position.
  • The one-dimensional interpolation procedure generalizes straightforwardly to any number of spatial dimensions.

Direct projection of the charge and current densities onto the grid

Charge and current densities are projected directly onto grid points by treating the projected quantity as the amount contained in the corresponding grid cell. The projection uses the recursive shape function and must be paired with matching interpolation shapes.

  • The projected quantity Q = (ρ, J) is assigned to a grid point by considering it as the charge or current contained in the surrounding cell.
  • Applying the recursive shape-function definition produces the direct projection expression.
  • Using the same shape function for interpolation and projection is mandatory to prevent unphysical self-force on quasi-particles.
Loading 1702.05128v1…