$π$-PIC: a framework for modular particle-in-cell developments and simulations

TL;DR

A Python-controlled framework designed to stimulate the dissemination and adoption of novel PIC developments by providing unified interfaces for accommodation, cross-testing, and comparison of PIC solvers and extensions written in Python or low-level languages like C++ and Fortran is presented.

Abstract

Recently proposed modifications of the standard particle-in-cell (PIC) method resolve long-standing limitations such as exact preservation of physically conserved quantities and unbiased ensemble down-sampling. Such advances pave the way for next-generation PIC codes capable of using lower resolution and fewer particles per cell, enabling interactive studies on personal computers and facilitating large-scale parameter scans on supercomputers. Here, we present a Python-controlled framework designed to stimulate the dissemination and adoption of novel PIC developments by providing unified interfaces for accommodation, cross-testing, and comparison of PIC solvers and extensions written in Python or low-level languages like C++ and Fortran. To demonstrate flexibility of proposed interfaces we present and test implementations of several PIC solvers, as well as of extensions that is capable of managing QED processes, moving-window, and tight focusing of laser pulses.

Figures (7)

• Figure 1: To the left, a simplified main script for simulation with $\pi$-PIC is shown, a full version can be seen in \ref{['app:example_script']}. Line 1 initiates the simulation instance. Line 2-5 adds particles and fields to the simulation. Line 6-7 defines functions which loops over fields and particles for saving or modification. All predefined decorator functions in the Python interface has a data_double and data_int as arguments. This allows the user to pass an address to a Numpy array which can be used to record the particle or field state, see line 10-11. Line 8 adds the extension to the advance call. During the advance call (line 9) the extensions and solver are applied to each cell and associated particles in a multi-threaded process. Line 9 advances the simulation state, executing the PIC solver and field solver and added extensions. The use_omp key word is used to activate multi-threading. The existence of multi-threaded processes is denoted as double arrow in the right scheme. Line 9 and 10 loops through the particles and field, passing an address to the double arrays particle_dd and field_dd, respectively, for saving the simulation state. The different instances of the code in the scheme to the right are color coded to match the Python commands that activates them.
• Figure 2: Panel (a) shows the initial field of the simulation: a Gaussian pulse is propagating at $45^\circ$ angle towards the boundary where the masking function $R(y)$ is shown in orange color scale. Panel (b) shows the difference in absorption function with different shape factors $\alpha$. The red curves corresponds to boundaries with maximal absorption as shown in panel (c). For these simulations the boundary size was set to four wavelengths. Panel (d) shows the dependence of the absorption on boundary size. The shape factor corresponding to $\Delta t \alpha=1$ was used for these simulations. The absorption was calculated as the decrease in total electromagnetic energy before and after interaction with the boundary.
• Figure 3: Electric field strength of the initial laser pulse after initialization within (a) $\Omega_B$ and (c) $\Omega_b$, along with the corresponding field states after a time step $R_0/c$ shown in (b) and (d), respectively. In (a), the shaded gray square indicates the extent of $\Omega_b$, while the surrounding unshaded squares represent equivalent coordinates within $\Omega_B$. Panel (e) shows the average runtime across $15$ runs for the initialization and propagation of the electromagnetic fields as a function of $R_0$.
• Figure 4: A two-stream instability by initializing opposite and equal currents flowing in a constant background field. The instability was simulated using two solvers: (left) electrostatic solver and (right) a electromagnetic solver (spectral field solver with Boris particle pusher). The density perturbation in phase space, $n$, is normalized to the initial particle density $n_0$. The electrostatic solver resolves the excitation of the two stream instability, but as the instability grows - the electrostatic simulation starts to gain total energy and eventually becomes invalid.
• Figure 5: Longitudinal electric field (a, c) and plasma density (b, d) for 1D LWFA simulations using $\pi$-PIC with EC and FB solver as well as Smilei with M4 solver. The black rectangles in (a) and (b) indicate the enlarged area shown in (c) and (d) respectively.