Runko: Modern multiphysics toolbox for plasma simulations

TL;DR

Runko is an open-source, modular plasma simulation framework in C++ and Python, enabling flexible, high-performance astrophysical plasma modeling across various computing platforms.

Contribution

It introduces a versatile, extendable multiphysics toolkit for plasma simulations with a hybrid design combining C++ efficiency and Python usability.

Findings

Successfully simulated relativistic Weibel instability in 1D

Modeled relativistic kinetic turbulence in 2D

Performed 3D simulations of collisionless shocks

Abstract

Runko is a new open-source plasma simulation framework implemented in C++ and Python. It is designed to function as an easy-to-extend general toolbox for simulating astrophysical plasmas with different theoretical and numerical models. Computationally intensive low-level kernels are written in modern C++ taking advantage of polymorphic classes, multiple inheritance, and template metaprogramming. High-level functionality is operated with Python scripts. The hybrid program design ensures good code performance together with ease of use. The framework has a modular object-oriented design that allows the user to easily add new numerical algorithms to the system. The code can be run on various computing platforms ranging from laptops (shared-memory systems) to massively parallel supercomputer architectures (distributed-memory systems). The framework supports heterogeneous multiphysics simulations in which different physical solvers can be combined and run simultaneously. Here we showcase the framework's relativistic particle-in-cell (PIC) module by presenting (i) 1D simulations of relativistic Weibel instability, (ii) 2D simulations of relativistic kinetic turbulence in a suddenly stirred magnetically-dominated pair plasma, and (iii) 3D simulations of collisionless shocks in an unmagnetized medium.

Tables

Table 1. Table 1: Time-advance loop of particle-in-cell module

Description	Routine name	Module/Solver	Eqs.	Comments
Grid syncronization
Communicate $\widehat{𝑩}$	mpi_b0	corgi		global communication
Communicate $\widehat{𝑬}$	mpi_e0	corgi		global communication
Update halo regions	upd_bc0	fields		local communication
First half of $\widehat{𝑩}$ advance
Advance $\widehat{𝑩}$	push_half_b1	fields/Propagator	(15)	$2$nd order FDTD†
Communicate $\widehat{𝑩}$	mpi_b1	corgi		global communication
Update $\widehat{𝑩}$ halo regions	upd_bc1	fields		local communication
Particle advance
Interpolate fields to particle locations	interp_em	pic/Interpolator		Linear scheme†
Advance particle $𝒖$ and $𝒙$	push	pic/Pusher	(17)–(18);	Boris scheme:† (72)-(78)
Second half of $\widehat{𝑩}$ advance
Advance $\widehat{𝑩}$	push_half_b2	fields/Propagator	(15)	$2$nd order FDTD†
Communicate $\widehat{𝑩}$	mpi_b2	corgi		global communication
Update $\widehat{𝑩}$ halo regions	upd_bc2	fields		local communication
$\widehat{𝑬}$ advance
Advance $\widehat{𝑬}$	push_e	fields/Propagator	(14)	$2$nd order FDTD†
Current $\widehat{𝑱}$ calculation
Calculate $\widehat{𝑱}$	comp_curr	pic/Depositer		ZigZag scheme†
Initialize $\widehat{𝑱}$ halo regions	clear_vir_cur	fields
Communicate $\widehat{𝑱}$	mpi_cur	corgi		global communication
Update $\widehat{𝑱}$ halo regions	cur_exchange	fields		local communication
Particle communications
Tag outflowing particles	check_outg_prtcls	pic
Pack them into a message	pack_outg_prtcls	pic
Communicate particles	mpi_prtcls	corgi		global communication
Unpack received messages	unpack_vir_prtcls	pic
Receive incoming particles	get_inc_prtcls	pic		local communication
Remove outflowing particles	del_trnsfrd_prtcls	pic
Clear global message arrays	del_vir_prtcls	pic
Deposit current $\widehat{𝑱}$ to $\widehat{𝑬}$
Add $\widehat{𝑱}$ to $\widehat{𝑬}$	add_cur	fields	(14)
I/O and analysis
Write output files	io	io

Table 2. Table 2: List of runko  modules and solvers

Module / Solver	Instance	Description	Eqs.	References
fields		Maxwell’s field equations module
xxxPropagator		Field solver
	FDTD2	$2$nd order FDTD solver	(14) and (15)	Yee (1966)
	FDTD4	$4$th order FDTD solver		Greenwood et al. (2004)
	FDTDGen	FDTD solver with free coefficients		Blinne et al. (2018)
xxxFilter		Current smoothing
	Binomial2	$3$-point (binomial) digital filter		Birdsall & Langdon (1985)
	Compensator2	$3$-point “compensator” filter		ibid.
pic		Particle-in-cell module
xxxPusher		Particle $𝒙$ and $𝒖$ propagation
	base class	Simple velocity Verlet propagator	(18)	Verlet (1967)
	BorisPusher	Relativistic Boris pusher	(72)-(78)	Boris (1970)
	VayPusher	Vay pusher		Vay (2008)
	HigueraCaryPusher	Higuera-Cary pusher		Higuera & Cary (2017)
	rGCAPusher	Guiding-center pusher		Bacchini et al. (2020)
xxxInterpolator		Field interpolator to particle location
	LinearInterpolator	Linear $1$st order interpolator
	QuadraticInterpolator	Quadratic $2$nd order interpolator
	CubicInterpolator	Cubic $3$rd order interpolator
	QuarticInterpolator	Quartic $4$th order interpolator
xxxDepositer		Current deposition
	ZigZag	First order ZigZag scheme	(79)-(83)	Umeda et al. (2003)
	ZigZag_2nd	$2$nd order ZigZag scheme		Umeda (2005)
	Esikerpov_2nd	$2$nd order Esikerpov scheme		Esirkepov (2001)
	Esikerpov_4th	$4$th order Esikerpov scheme		ibid.

Equations

$${\gamma }^2=1+(u/c{)}^2=(1-(v/c{)}^2{)}^{-1},$$

$$\frac{\partial f_{\mathrm{s}}}{\partial t}+v\cdot {\nabla }_x{f}_{\mathrm{s}}+a_{\mathrm{s}}\cdot {\nabla }_u{f}_{\mathrm{s}}=\mathcal{C},$$

$$a_{\mathrm{s}}\equiv {\mathrm{d}}_{t}u=\frac{q_{\mathrm{s}}}{m_{\mathrm{s}}}\left(E+\frac{v}{c}\times B\right)=\frac{q_{\mathrm{s}}}{m_{\mathrm{s}}}\left(E+\frac{u}{\gamma c}\times B\right),$$

$${\rho }_{\mathrm{s}}\equiv q_{\mathrm{s}}\int f_{\mathrm{s}}\mathrm{d}u.$$

$$J_{\mathrm{s}}\equiv q_{\mathrm{s}}\int f_{\mathrm{s}}v\mathrm{d}u=q_{\mathrm{s}}\int f_{\mathrm{s}}\frac{u}{\gamma }\text{leavevmode}\mathrm{d}u.$$

$$\nabla \cdot E=4\pi \rho ,$$

$$\nabla \cdot B=0,$$

$$\nabla \times E=-\frac{1}{c}\frac{\partial B}{\partial t},$$

$$\nabla \times B=\frac{4\pi }{c}J+\frac{1}{c}\frac{\partial E}{\partial t}.$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot J=0.$$

$$\hat{c}=c\frac{\Delta t}{\Delta x},$$

$$\hat{E}$$

$$\hat{B}$$

$$\Delta [\hat{E}{]}_t=\hat{c}\Delta [\hat{B}{]}_x-\hat{J}$$

$$\Delta [\hat{B}{]}_t=-\hat{c}\Delta [\hat{E}{]}_x,$$

$$m\frac{\mathrm{d}u}{\mathrm{d}t}=q\left(E+\frac{v}{c}\times B\right).$$

$$\Delta [\hat{u}{]}_t=\frac{\hat{q}}{\hat{m}}\left(\hat{E}+\frac{\hat{v}}{\hat{c}}\times \hat{B}\right).$$

$$\Delta [x{]}_t=\frac{u\hat{c}}{\gamma }.$$

$$\frac{\delta }{{\omega }_{\mathrm{p}}}={\beta }_{\mathrm{b}}\sqrt{\frac{2}{{\gamma }_{\mathrm{b}}}},$$

$$\sigma =\frac{B_0^2}{4\pi n_0\bar{\gamma }{m}_{\mathrm{e}}{c}^2},$$

$$\delta B_x$$

$$\delta B_y$$

$${\beta }_{mn}=2\frac{\delta B}{N\sqrt{m^2+n^2}},$$

$$\nabla \cdot E$$

$$\frac{\partial E}{\partial t}$$

$$\nabla \cdot B$$

$$\frac{\partial B}{\partial t}$$

$$\nabla \cdot E$$

$$\frac{\partial E}{\partial t}$$

$$\nabla \cdot B$$

$$\frac{\partial B}{\partial t}$$

$$\frac{\mathrm{d}u}{\mathrm{d}t}=\frac{q}{m}\left(E+\frac{v}{c}\times B\right).$$

$$\Delta [\hat{E}{]}_t$$

$$\Delta [\hat{B}{]}_t$$

$$\Delta [\hat{u}{]}_t$$

Full text

1 Introduction

Ever since the introduction of computers, numerical simulations have been used to study the nonlinear behavior of plasma (see, e.g., Buneman 1959; Dawson 1962, 1964). In the early days, the research was mainly motivated by studies of basic plasma instabilities and confinement in fusion experiments, but the method of computationally solving the dynamics of charged particles quickly gained popularity also in understanding plasma in space (see, e.g., Tanaka 1983; Langdon et al. 1988; Buneman et al. 1992).

An important factor in accelerating the use of plasma simulations is the ever-increasing computational speed and number of processors. We have recently arrived in the petaflop ($10^{15}$ floating-point operations per second, FLOPS) supercomputing era, and have started to pave our way toward exascale computing systems ($10^{18}$ FLOPS). This technological surge (and the ingenuity of the researchers themselves) has enabled us to shed light on many longstanding issues in high-energy astrophysics, including shocks (e.g., Frederiksen et al. 2004; Spitkovsky 2005), reconnection (e.g., Zenitani & Hoshino 2001; Cerutti et al. 2012; Kagan et al. 2013; Sironi & Spitkovsky 2014; Werner et al. 2015), pulsar magnetospheres (e.g., Philippov & Spitkovsky 2014; Chen & Beloborodov 2014; Cerutti et al. 2015), kinetic turbulence (e.g., Howes et al. 2008; Servidio et al. 2014; Zhdankin et al. 2017b; Comisso & Sironi 2018; Nättilä & Beloborodov 2021), and nonlinear wave interactions (e.g., Nättilä & Beloborodov 2022).

In order to keep riding this technological wave, successful numerical codes need to be highly scalable and modern in order to take advantage of the newest supercomputing architectures. These kinds of exascale-ready features include avoidance of global communication patterns and support of heterogeneous computing nodes with varying resources, such as a different number of cores/accelerators or varying amount of memory per core (see also Nordlund et al. 2018, for a discussion). The presented runko framework was developed from scratch to support these features.

An additional feature of the plasma simulations is that there are various theoretical frameworks for modeling the dynamics of the system, whereas the underlying numerical algorithms all rely on techniques based on computational particles and multi-dimensional grids. The fully kinetic description is based on sampling the distribution function with superparticles and representing the electromagnetic fields on a grid, known as the particle-in-cell (PIC) approach (Hockney & Eastwood 1981; Birdsall & Langdon 1985).111Alternatively, such a system can be modeled by solving the partial differential equations directly on a grid (known as the Vlasov method; Cheng & Knorr 1976).

In the case of a strong (background) magnetic field, it is possible to average over the gyrorotation of the particle’s gyro-orbit and use the gyrokinetic description where, instead of particles, computational “rings” are considered (see, e.g., Howes et al. 2006). Alternatively, the kinetic systems can be significantly simplified by considering only the moments of the distribution function; this macro-scale description is related to the magnetohydrodynamic (MHD) models (see, e.g., Brandenburg 2003). The MHD description itself can be divided into multiple categories with different closure approximations (e.g., ideal, visco-resistive) or the evolution of multiple fluid species (known as multifluid description). There also exists a class of hybrid methods that combine both kinetic and magnetohydrodynamic models, for example by treating electrons with fluid description and protons with the kinetic description (see, e.g., Winske et al. 2003; Valentini et al. 2007; Haggerty & Caprioli 2019). From the numerical perspective, all of these methods, in the end, only require the efficient manipulation of particles and/or physical quantities that are stored on grids.

It is possible to take advantage of this structural similarity of the various formulations with well-designed software abstractions. This is the second core idea of the runko simulation framework, which in essence is just a collection of data structures for storing, manipulating, and communicating particles and grids efficiently. The framework relies on modern computer science methods and paradigms that provide a sufficient degree of problem generalization to technically enable easy introduction of physics knowledge into the simulations. Technically this is achieved by polymorphic multiple inheritance as supported by the use of modern C++. The concept of object inheritance also nicely reflects the mathematical way of constructing hierarchical models that build on the previous one. Another level of abstraction is obtained by using template metaprogramming that allows, for example, the parameterization of the dimensionality of the system, again closely reflecting the mathematical form of the problem. Finally, in order to balance the computational efficiency and the user-friendliness of the code, the fast low-level C++ classes are exposed to high-level Python3 scripting language.

The complete framework is provided as open-source software for the community. It is available from a GitHub repository222https://github.com/natj/runko 

and the code design and structure are presented in detail in this article. Here we also present the PIC code of the framework and describe the implementation of the solvers needed to operate it. First, the theoretical background of the kinetic modeling of plasma is briefly reviewed in Sect. 2. The numerical implementation of the framework is then discussed in Sect. 3. As an example, the use of the new PIC code is demonstrated on astrophysical plasma problems. A discussion of new features and the future directions of the framework are given in Sect. 5. Finally, the content is summarized in Sect. 6.

2 Theory

2.1 Kinetic plasma theory

The kinetic plasma theory we adopt here is based on the special relativistic formulation of the Vlasov–Boltzmann equation. The spatial coordinate location vector is given by $\bm{\mathit{x}}\equiv(x,y,z)$ and coordinate time is measured by $t$. Coordinate velocity (three-velocity) is given by $\bm{\mathit{v}}\equiv d_{t}\bm{\mathit{x}}$ and the individual Cartesian components of the velocity are denoted as $\bm{\mathit{v}} =(v_{x},v_{y},v_{z})$. Proper time $\tau$ (measured with a co-moving clock) is connected to the coordinate time via the Lorentz factor $\gamma\equiv d_{\tau} t$. The proper velocity (spatial components of the four-velocity) is $\bm{\mathit{u}} \equiv d_{\tau}\bm{\mathit{x}}=\gamma\bm{\mathit{v}}$. The Lorentz factor and the velocities are connected by the expression

[TABLE]

where $c$ is the speed of light, $u=|\bm{\mathit{u}}|$, and $v=|\bm{\mathit{v}}|$. Acceleration is denoted with $\bm{\mathit{a}} \equiv d_{t}\bm{\mathit{u}}$.

A six-dimensional phase-space density distribution for particle species s is given by $f_{\mathrm{{s}}}\equiv f_{\mathrm{{s}}}(\bm{\mathit{x}},\bm{\mathit{u}};t)$. Thus, $f_{\mathrm{{s}}}\,d^{3}x\,d^{3}u$ is the number of particles in the six-dimensional differential phase space volume between $\bm{\mathit{x}}$, $\bm{\mathit{u}}$ and $\bm{\mathit{x}}+d\bm{\mathit{x}}$, $\bm{\mathit{u}}+d\bm{\mathit{u}}$.

The evolution of $f_{\mathrm{{s}}}$ is governed by the relativistic Boltzmann–Vlasov equation

[TABLE]

where $\nabla_{\bm{\mathit{x}}}=\frac{\partial}{\partial\bm{\mathit{x}}}$ and $\nabla_{\bm{\mathit{u}}}=\frac{\partial}{\partial\bm{\mathit{u}}}$ are the spatial and momentum parts of the differential operator $\nabla$, respectively. The term on the right-hand side, defined as $\mathcal{C}\equiv\partial_{t}f_{\mathrm{{s}}}\leavevmode\nobreak\ |_{\mathrm{coll}}$, is the collision operator. For the Vlasov equation $\mathcal{C}=0$ (i.e., the plasma is collisionless). Acceleration of a charged particle $\bm{\mathit{a}}_{\mathrm{s}}$ is governed by the Lorentz force

[TABLE]

where $\bm{\mathit{E}}$ and $\bm{\mathit{B}}$ are the electric and magnetic fields, $q_{\mathrm{s}}$ is the charge, and $m_{\mathrm{s}}$ is the mass of the particle of species s.

The moments of the distribution function define macroscopic (bulk) quantities of the plasma. The zeroth moment of the distribution function $f_{\mathrm{{s}}}$ defines the charge density of species s as

[TABLE]

The total charge density is $\rho=\sum_{\mathrm{s}}\rho_{\mathrm{s}}$. The first moment defines the current (charge flux) vector as

[TABLE]

The total current is $\bm{\mathit{J}}=\sum_{\mathrm{s}}\bm{\mathit{J}}_{\mathrm{s}}$.

2.2 Maxwell’s equations

Evolution of electric field $\bm{\mathit{E}}$ and magnetic field $\bm{\mathit{B}}$ is governed by Maxwell’s equations. These are Gauss’s law

[TABLE]

Gauss’s law for magnetism

[TABLE]

Ampere’s law

[TABLE]

and Faraday’s law

[TABLE]

Charge conservation follows from these equations by taking a divergence of Eq. (9) and substituting Eq. (6) to get

[TABLE]

2.3 Dimensionless equations

We now describe the actual dimensionless equations that are solved numerically. A unit system similar to that used in Tristan and Tristan-MP(Buneman et al. 1993; Spitkovsky 2005; Sironi et al. 2013) is used. The derivation and more thorough discussion of this unit system is given in Appendix A.

Many quantities, such as the electromagnetic fields, are defined on a lattice (mesh). The discrete location of a point in this case is given as $\bm{\mathit{x}}^{(i,j,k)}\equiv\bm{\mathit{x}}(i,j,k)=(i\Delta x,j\Delta y,k\Delta z),$ where $i,j,k$ are the grid indices and $\Delta x$, $\Delta y$, and $\Delta z$ are the grid point distances in each Cartesian dimension. Similarly, discretized time is given as $t^{n}\equiv t(n)=n\Delta t$. Cells of the lattice are taken to have a fixed cube geometry: $\Delta x=\Delta y=\Delta z$. The Courant (or Courant-Friedrichs-Lewy; Courant et al. 1928) number is defined as

[TABLE]

which represents a maximum numerical velocity of a signal traveling on the grid. For explicit time integration schemes $\hat{c}\leq 1$.

Electromagnetic fields are normalized with $B_{0}$ as $\bm{\mathit{\hat{E}}}=\bm{\mathit{E}}/B_{0}$ and $\bm{\mathit{\hat{B}}}=\bm{\mathit{B}}/B_{0}$. Similarly, current density is normalized with $J_{0}$ as $\bm{\mathit{\hat{J}}}=\bm{\mathit{J}}/J_{0}$. The values of these normalization factors are selected such that the numerically solved equations appear as $\Delta x=\Delta t=1$. This also means that the grid indices $i$, $j$, and $k$ have the same numerical value as location $\bm{\mathit{x}}$.

2.3.1 Electromagnetic field module

The time evolution of electromagnetic fields is handled with a finite-difference time-domain (FDTD) method. A Yee lattice (Yee 1966) is used for the $\bm{\mathit{E}}$ and $\bm{\mathit{B}}$ fields such that they are staggered both in space and in time

[TABLE]

where $\bm{\mathit{\hat{E}}}$ is located in the middle of the cell sides and $\bm{\mathit{\hat{B}}}$ in the center of the cell faces. This makes it easy to calculate the curl of the fields in the subsequent equations.

In the time domain we update the $\bm{\mathit{\hat{E}}}$ and $\bm{\mathit{\hat{B}}}$ fields with discrete forms of Eqs. (8) and (9) given as

[TABLE]

and

[TABLE]

where $\Delta[Q]_{t,\bm{\mathit{x}}}$ is the time differentiation or curl of a variable $Q$ without the $\Delta x$ or $\Delta t$ denominator. The only normalization factor entering these equations is the Courant velocity $\hat{c}$ since everything else is absorbed in $B_{0}$ and $J_{0}$. There is no need to solve Eqs. (6) and (7) if the charge-conserving scheme is used together with the Yee staggering of the fields (see Appendix B for details).

2.3.2 Particle-in-cell module

The main purpose of the pic module is to update particle velocities according to the Lorentz force. We express the spatial component of the four-velocity $u$ in units of $c$. The particle species specifier s is omitted in this section for brevity.

The Lorentz force acting on a charged particle is

[TABLE]

We simplify this expression again to make it appear unitless. We express the charge and mass as $q=q_{0}\hat{q}$ and $m=m_{0}\hat{m}$. The right-hand side of the Lorentz force equation is then

[TABLE]

The actual velocity update is done in parts since $\bm{\mathit{\hat{E}}}$ and $\bm{\mathit{\hat{B}}}$ are staggered in time. As an example of particle momentum update, the relativistic Boris scheme (Boris 1970) is presented in detail in Appendix C. In addition to the particle velocity update routines, we need a scheme to interpolate the (staggered) $\bm{\mathit{\hat{E}}}$ and $\bm{\mathit{\hat{B}}}$ fields to the arbitrary particle locations. This interpolation is typically implemented using a linear volume-averaging interpolation scheme although higher-order schemes can also be utilized.

After the velocity update, we advance the particle’s location. The particle’s position $\bm{\mathit{x}}$ is advanced in time with

[TABLE]

The described particle propagation scheme is second order in time. The velocity is stored in the middle of the time steps and updated as $\bm{\mathit{u}}_{n+\frac{1}{2}}=\bm{\mathit{u}}_{n-\frac{1}{2}}+\Delta[\bm{\mathit{\bm{\mathit{\hat{u}}}}}]_{t,\,n}$; the location is stored at integer time steps and updated as $\bm{\mathit{x}}_{n+1}=\bm{\mathit{x}}_{n}+\Delta[\bm{\mathit{x}}]_{t,\,n+\frac{1}{2}}$.

Finally, when a charged particle moves on the grid it creates a current $\bm{\mathit{\hat{J}}}$ that induces an electric field via Eq. (9). We use the charge-conserving ZigZag current deposition scheme to deposit the current from the computational particles into the mesh (Umeda et al. 2003). This scheme is summarized in Appendix D.

3 Numerical implementation

In this section we now discuss the general design principles and structure of the framework itself. runko is a hybrid C++/Python code: low-level computationally intensive kernels are written in C++, whereas the high-level functionality is operated by Python scripts. The design heavily relies on an object-oriented programming paradigm (supported naturally by both C++ and Python), where data containers (attributes) and functions (methods) are grouped together into objects that interact with each other. In runko we call a group of these objects modules (Sect. 3.2) and the objects themselves are called solvers (Sect. 3.3). The solvers are operated by different kinds of drivers (Sect. 3.4). This kind of hybrid use of C++ and Python ensures a fast and well-optimized low-level code, whereas the Python scripts allow for ease of use and straightforward extensibility of the code.

The low-level C++ kernels rely heavily on template metaprogramming features of modern C++ (conforming currently to C++$14$ standard333The C++$14$ and $17$ standards are generally considered as “transition standards” bridging the gap between the older C++$11$ and the newer C++$20$ version. In the future, runko will transition to fully conform to the C++$20$ standard as soon as sufficient HPC compiler support can be guaranteed. At the moment, many of the newer C++$20$ template metaprogramming features are provided together with runko itself, implemented using C++$11$. ). In template metaprogramming a code template is provided to the compiler that then generates the actual source code based on it, given some extra information on how to fill the template parameters.444We use the standard C++ syntax to present template classes: an object A with a template parameter X is given as A<X>.

Many of the implemented C++ kernels take advantage of this by only having a $D$-dimensional general template of the algorithm. This template of an algorithm is then specialized to the required spatial dimension (typically $D=1$, $2$, or $3$) by the compiler during compile time. This translates to a more bug-free code as typically only one implementation of the algorithm is required.

The low-level C++ kernels are operated by user-friendly Python scripts. Technically, the C++ kernels are exposed to Python using the Pybind11 library (Jakob et al. 2017). All the exposed C++ objects behave like native Python objects. No unnecessary memory copies are performed because Pybind11 manages memory pointers. After installing runko these objects are available to Python by importing them from the pyrunko package (and pycorgi, as described in Sect. 3.1). Not only can these classes be easily used from the Python scripts, but they also support inheritance in the case the user wants to extend and modify the objects. This allows for fast prototyping as new classes can first be developed in Python before implementing them in C++ (if better performance is needed). These Python scripts also leverage Numpy (van der Walt et al. 2011) and Scipy for performance and convenience, together with Matplotlib (Hunter 2007) for visualization.

3.1 Grid infrastructure

runko is built on top of the massively parallel corgi555https://github.com/natj/corgi C++ template library (Nättilä, in prep.). corgi is a modern parallelized grid infrastructure library that provides the (spatial) simulation grid. It is designed to be run on distributed memory systems (i.e., computing clusters and supercomputers) that require explicit inter-node communication. All inter-node communications are handled by the corgi library using the Message Passing Interface (MPI) library.666https://github.com/natj/mpi4cpp.

Hereafter, one MPI process is called a rank, to make a distinction with a processor because a rank does not necessarily have to coincide with a processor.

corgi uses a patch-based super-decomposition strategy to partition the simulation grid between different ranks. This means that the full simulation grid is split into small (preferably continuous) subregions called tiles. One tile can, for example, host a $10\times 10$ subgrid in a 2D simulation. If the complete simulation domain is composed of a $10^{4}\times 10^{4}$ lattice, it can be tessellated with $10^{3}\times 10^{3}$ of these tiles.

The implementation of a physical simulation code works by inheriting all the grid infrastructure methods from a corgi::Tile<D> template class. This derived template class then automatically holds all the relevant physical attributes and methods needed to perform the simulation communications. The tile object provided by corgi is a template class that is parameterized with the spatial dimension $D$. Specializing this dimension parameter to any given value sets the simulation grid dimensionality.

All the tiles are stored in a container provided by corgi called a grid. Each program rank has its own grid object corgi::Grid<D> (specialized for $D$) that holds a varying number of tiles. This allows us to decompose the grid so that each rank will only hold and process a small part of the simulation domain (i.e., set of tiles). The grid class is responsible for providing all the spatial neighborhood information of the tiles and executing the inter-tile communications. These include both local shared-memory intra-rank communications and global inter-rank MPI communications. In practice, these methods are typically used to keep the halo regions of the tiles up to date.

There are no restrictions on the shape of the tile boundaries between different ranks. This allows more complex tile ownership to be used between ranks, which aims to maximize the volume and minimize the surface area, like a honeycomb tessellation in 2D. This configuration then translates directly to less communication needed between ranks.

Both the corgi::Tile and corgi::Grid C++ classes are exposed to Python. Dimensionality (i.e., the $D$ template parameter) is automatically set by loading a correct subpackage from the runko (and corgi) Python libraries: options are oneD for one-dimensional simulations ($N_{x}$ sized tiles; $\mathcal{N}_{x}$ dimensional global grid configuration), twoD for two-dimensional simulations ($N_{x}\times N_{y}$ sized tiles; $\mathcal{N}_{x}\times\mathcal{N}_{y}$ dimensional global grid configuration), and threeD for three-dimensional ($N_{x}\times N_{y}\times N_{z}$ sized tiles; $\mathcal{N}_{x}\times\mathcal{N}_{y}\times\mathcal{N}_{z}$ dimensional global grid configuration). Here $N_{x,y,z}$ is the subgrid resolution of tiles and $\mathcal{N}_{x,y,z}$ is the simulation grid resolution in units of the tiles.

3.2 Modules

The actual physics of the simulation is implemented in derived tiles that are based on the corgi base tiles. We call these tiles and all the related objects and functions modules.777Technically, a module corresponds to a namespace in C++ and a subpackage in Python.

An important part of the framework design is that different modules can be combined with polymorphic multiple inheritance (i.e., deriving methods and attributes from multiple other modules). This allows an efficient re-use of different modules when designing a new module with similar or overlapping physics concepts.

As an example, the PIC code (Sect. 2.3.2) uses corgi::Tile as its first base class to inherit the grid infrastructure and fields::Tile as its second base class to inherit the electromagnetic FDTD routines (Sect. 2.3.1). The pic module itself only defines the particle containers and methods to manipulate particle movement. The two modules are connected by the current $\bm{\mathit{\hat{J}}}$ that is created by the moving charged particles. This current then induces a change in the $\bm{\mathit{\hat{E}}}$ and $\bm{\mathit{\hat{B}}}$ fields that affects the Lorentz force experienced by the particles (see Fig. 1).

3.3 Solvers

Each module can contain several solvers that operate on and/or change the internal attributes of the module’s tiles. Solvers are created by first defining an interface with an abstract class. The actual solver implementations that can be instantiated (i.e., the concrete classes) can then be implemented by creating a derived class from this abstract interface. These abstract interface classes are bound to Python objects via the “trampoline” classes. This automates the binding of any user-defined concrete solver implementations. A typical requirement of any solver is that it should provide a solve method that takes (at least) the module’s tile as an input.

As an example, the pic module defines an abstract Interpolator solver interface. The purpose of this solver object is to interpolate the $\bm{\mathit{\hat{E}}}$ and $\bm{\mathit{\hat{B}}}$ fields, defined on a lattice, to arbitrary particle locations. One of the currently provided Interpolator instances is the Linear solver that implements a linear first-order interpolation scheme to approximate the field values between the grid points. Thanks to this design pattern, it is very easy to extend the module to include other interpolation methods, such as a second-order polynomial or a spectral interpolation. The framework is agnostic about how the job of the solver is actually fulfilled; it is enough that it conforms to the abstract interface class standard, after which it can be automatically incorporated into the module routines.

3.4 Drivers

The actual simulation codes in runko are executed by the drivers. Drivers are short Python scripts that use the C++ bindings of the tiles and solvers. They typically consist of an initialization section in the beginning where tiles of some specific module are loaded into the corgi grid container. After this, the data attributes of the tiles are initialized. The drivers also include the main time propagation loop where solvers are applied to the tiles in the correct order. As an example, Table 1 presents the main time integration loop of the PIC code.

The main advantage of using the low-level C++ classes via the Python driver scripts is the ease of use. For example, various complex initialization routines are rendered much easier to implement in a high-level language that also supports data visualization. It also opens up the possibility of performing the analysis (or at least part of it) on the fly during the simulation.

3.5 Currently implemented modules

The currently implemented modules include the fields and pic modules. These and the corresponding solvers are listed in Table. 2.

The fields module (Sect. 2.3.1) is responsible for advancing the electromagnetic fields via the Maxwell’s equations. This module inherits its spatial dimensionality $D$ from the corgi tiles. Internally it has a mesh container that stores the field components as $D$-dimensional lattices. It also holds methods for advancing these fields in time using Eqs. (14) and (15). The module is agnostic about the form of the incoming current; it can be coupled to any other module to provide this closure.

The pic module (Sect. 2.3.2) handles the particle propagation and provides the current closure for the fields module (see also Fig. 1). Attributes of the electromagnetic field lattices are inherited from the fields::Tile. In addition to the spatial dimensionality template parameter $D$, the pic module has another template parameter for the velocity space dimensionality $V$ (i.e., pic::Tile<D,V>). By specializing the $V$ template parameter, the solvers acting on the pic::Tile take into account the corresponding velocity space dimensionality. For example, setting $V=2$ equals manipulating the velocity vector as $\bm{\mathit{u}}=(u_{x},u_{y})$, whereas $V=3$ modifies a full three-dimensional velocity vector $\bm{\mathit{u}}=(u_{x},u_{y},u_{z})$. The pic module also defines various solvers that operate on the pic::Tiles (see Table 2 for all the solvers implemented in the framework).

The final important module in the current framework is the io module, which is responsible for the data input and output functionality. It is currently more detached from the physical implementations providing only solvers that act on the tiles. These solvers can be categorized into two distinct types of input readers and output writers. The input readers load data into the tiles from storage files, whereas the output writers write data from the tiles to the disk. Currently implemented readers and writers are capable of operating on full simulation snapshot restart files using the hierarchical data format (HDF5). In addition, both field and particle classes have more storage-efficient writers that reprocess the data on the fly during the simulation and only write a smaller subsample of the full data into disk.

4 Problem showcase

In this section we showcase the new PIC code implemented into the runko framework. First, we demonstrate the code validity by comparing the simulation results of the relativistic Weibel instability to the known analytical growth rates. Next, we demonstrate the same instability on a larger scale by presenting results from an unmagnetized collisionless shock simulation. Finally, we measure the numerical performance of the code by using results from relativistic turbulence simulations.

4.1 Relativistic Weibel instability

We start showcasing the code capabilities by first simulating a standard astrophysics-motivated plasma instability, commonly known as the Weibel instability (Weibel 1959; Fried 1959; Bret et al. 2010). The relativistic Weibel instability in a collisionless plasma is triggered when two counter-streaming plasma populations start to experience small charge-dependent deflections, seeded by the noise. The resulting modes are electromagnetic so perturbations in both $\bm{\mathit{E}}$ and $\bm{\mathit{B}}$ are important. Therefore, the presented simulations also help us probe the validity of the solutions to the full Vlasov–Maxwell equations.

4.1.1 Problem setup

We simulate the linear growth phase of the instability by following two counter-streaming pair plasma populations in a 2D simulation domain. The plasma is modeled with $32$ particles per cell per species, and is initialized to have a small thermal spread of $\theta=kT/m_{\mathrm{e}}c^{2}=10^{-5}$ and to flow relativistically along $\hat{\bm{\mathit{x}}}$ as $\bm{\mathit{u}}=\pm\gamma_{\mathrm{b}}\hat{\bm{\mathit{x}}}$. Electromagnetic fields are propagated in time with the standard FDTD2 solver, current is deposited with the ZigZag scheme, particle interpolation is done with the LinearInterpolator method, and particle propagation is performed with the BorisPusher solver. The simulation domain is captured with a grid of $320\times 80$ cells and the skin depth is resolved with $10$ cells. We use a standard time step of $\Delta t=0.45\Delta x$. We use symmetric beams with Lorentz factors of $\gamma_{\mathrm{b}}=3$, $10$, $30$, and $100$.

4.1.2 Comparison to analytical estimates

The analytic linear growth rates of the Weibel instability are known for symmetric beams. Bret et al. (2010) presents the filamentation instability growth rate $\delta$ as

[TABLE]

where $\gamma_{\mathrm{b}}=(1-\beta_{\mathrm{b}}^{2})^{-1/2}$ is the bulk Lorentz factor of the beams, $\omega_{\mathrm{p}}=\sqrt{4\pi n_{\mathrm{e}}e^{2}/\gamma_{\mathrm{b}}m_{\mathrm{e}}}$ is the plasma frequency, and $n_{\mathrm{e}}$ is the plasma number density (of electrons). We compare the simulation results to these theoretical estimates in Fig. 2 and find the linear instability growth rate to be perfectly matched with the expectations.

These simulations are not particularly demanding and can therefore serve as a helpful introduction to using the code. The discussed setup is also general enough that it could be easily reformulated to study other analytically intractable plasma beam instabilities.

4.2 Formation of collisionless shocks

The relativistic Weibel filamentation setup can be expanded to a realistic 3D simulation and used to model the plasma physics of collisionless shocks. The shock is formed by the growth and nonlinear saturation of the same instability.

We find this setup to be one of the most demanding numerical tests for PIC codes because the numerical algorithms need to be able to both sustain the relativistic upstream bulk flow and correctly model the kinetic shock physics. The upstream is unstable to the numerical Cherenkov instability (Godfrey 1974) and is, in addition, sensitive to any numerical error in the current deposition; even tiny perturbations in the relativistically advected upstream plasma can trigger an isotropization of the bulk flow. The shock front, on the other hand, is sustained by a strong electric current that is sensitive to charge conservation. We use this setup to stress test the full 3D PIC routines and their correctness.

4.2.1 Problem setup

The numerical setup consists of a elongated shock “tube” along the $\hat{\bm{\mathit{x}}}$-axis and a reflecting wall on the left side of the domain. Initially, we set all the plasma in the domain to move from right to left with a velocity $\bm{\mathit{u}}=-\Gamma\beta\hat{\bm{\mathit{x}}}$. The plasma is reflected from the left simulation domain and forms a counter-propagating beam. The beam is unstable to the Weibel instability that will quickly stop the right-moving component. This causes the formation of a collisionless shock.

We simulate the shock in an elongated box with a size $L_{\parallel}\times L_{\perp}^{2}=512\times 64^{2}(c/\omega_{\mathrm{p}})^{3}$ using ten cells per skin depth so that the full domain is covered with a grid of $5120\times 640^{2}$. The plasma flows with a relativistic velocity $\Gamma=10$ and is modeled with 64 particles per cell per species. Electromagnetic fields are propagated in time with the fourth-order FDTD4 scheme. The resulting current is calculated using the ZigZag scheme, particle interpolation is done via the LinearInterpolator, and the current is smoothed with eight passes of the binomial filter (Binomial2). We use a time step of $\Delta t=0.45\Delta x$. The pair plasma has an initial thermal spread of $\theta=kT/m_{\mathrm{e}}c^{2}=2\times 10^{-5}$ and is updated using the HigueraCaryPusher. We start the simulation with no background magnetic field, which corresponds to an unmagnetized shock.

The plasma is injected into the domain from a receding right wall that moves at the speed of light, from left to right; this enables longer simulations both because the particles are not “alive” as long as they normally would be (suppressing numerical Cherenkov instability) and by saving computational time (since we do not need to initially propagate all the particles in the domain). The shock front is adaptively tracked on the fly, and only the high-resolution grid data around it is saved. The grid is decomposed into cyclic MPI domains of length $1280$ cells in the $\hat{\bm{\mathit{x}}}$ direction (i.e. five cycles per the full box); this results in a more even memory load for the simulation.

4.2.2 Physical interpretation and numerical tests

A collisionless shock is quickly formed when the simulation is started. The shock is mediated by a formation of Weibel filaments that grow, saturate, and merge; this stops the counter-streaming plasma (see Fig. 3).

The simulated shock structure and dynamics are well described by the theoretical expectations and previous numerical experiments (e.g., Sironi & Spitkovsky 2009; Plotnikov et al. 2018). The filaments form and induce the shock in about $100\omega_{\mathrm{p}}^{-1}$. The upstream flow is compressed into a denser downstream with a compression ratio $n_{\mathrm{down}}/n_{\mathrm{up}}\approx 4.1\pm 0.15$, and the shock propagates into the upstream with a mildly relativistic velocity, $\beta_{\mathrm{sh}}\approx 0.32\pm 0.03$ (measured for $200\leq t\omega_{\mathrm{p}}\leq 1000$). Unmagnetized MHD jump conditions, in comparison, predict, $n_{\mathrm{down}}/n_{\mathrm{up}}=1+(\Gamma+1)/[\Gamma(\gamma_{\mathrm{AD}}-1])=4.3$ for the compression ratio, and $\beta_{\mathrm{sh}}=(\gamma_{\mathrm{AD}}-1)(\Gamma-1)/\Gamma)=0.3$, given the 3D adiabatic index $\gamma_{\mathrm{AD}}=4/3$ (see, e.g., Plotnikov et al. 2018). The shock ramp itself has a length of roughly $\approx 50c/\omega_{\mathrm{p}}$ and is accommodated by an increase in the electric field energy density just ahead of the moving front.

We also specifically tested that the upstream remains smooth and stable for an extended period of time $\gtrsim 1000\omega_{\mathrm{p}}^{-1}$ with electric and magnetic field fluctuations of $\delta E$, $\delta B\lesssim 10^{-6}$ (in absolute code units; the value corresponds to floating-point accuracy). This is a stringent test that can be used to verify that the particle-to-grid and grid-to-particle routines are all correctly implemented, and no charge or electric current is lost; in the opposite case, any numerical error will quickly grow and isotropize the flow. Similarly, the formation of the shock front supports the notion that all the 3D routines are correct. We note that the shock simulation is a particularly demanding code test because of the above details.

4.3 Decaying relativistic kinetic turbulence

In this section we use a relativistic decaying turbulence setup to showcase the code performance and scalability metrics. The setup is presented and discussed in detail in Comisso & Sironi (2018, 2019); Nättilä & Beloborodov (2021). Similar driven turbulence simulations are performed in Zhdankin et al. 2017a, b, among others.

4.3.1 Problem setup

We consider a uniform pair plasma that is suddenly stirred. Due to this large-scale driving motion, a turbulent cascade quickly develops. We study the 2D3V setups, where the cascade forms in the planar $x-y$ direction and is captured by the spatial 2D simulation domain. All three velocity components are considered, however, including the out-of-the-plane $z$ direction. To physically mimic confinement of the cascade to a 2D plane, a guide field is imposed in the $z$ direction as $\bm{\mathit{B}}_{0}=B_{0}\hat{\bm{\mathit{z}}}$ ($\hat{\bm{\mathit{z}}}$ is the unit vector along the $z$-axis).

Presence of the strong guide-field renders the plasma magnetically dominated. The plasma magnetization parameter we use is

[TABLE]

where $n_{0}$ is the initial number density and $\bar{\gamma}$ is the mean thermal Lorentz factor of the plasma. Here we focus only on a strongly magnetized case of $\sigma=10$. Initially, the plasma is defined to be warm with a temperature of $\theta=kT/m_{\mathrm{e}}c^{2}=0.3$ ($\bar{\gamma}\approx 1.6$).

The spatial resolution of our fiducial simulation is $5120^{2}$ cells ($512\times 512$ tiles with internal tile resolution of $10\times 10$). The initial plasma skin depth $c/\omega_{\mathrm{p}}$ is resolved with five cells. This determines our box size as $L\approx 1000c/\omega_{\mathrm{p}}$. We use $64$ particles per cell per species, totaling $128$ particles per cell for both electron and positron species. In total, we then propagate $\sim\!10^{9}$ particles in the simulation.

The sudden large-scale stirring is mimicked by perturbing the planar magnetic field $\bm{\mathit{B}}_{\perp}$ with strong $\langle B_{\perp}^{2}\rangle\equiv\delta B^{2}\sim B_{0}^{2}$ uncorrelated fluctuations in the $x$ and $y$ directions following Comisso & Sironi (2018). This drives strong turbulence with semi-relativistic bulk velocities. The forcing modes, $m,n\in\{1,\ldots,8\}$, are initialized with equal power as

[TABLE]

where $k_{m}=2\pi m/L$ and $k_{n}=2\pi n/L$ are the wave numbers along the $x$ and $y$ dimensions, respectively. The forcing is uncorrelated so $\phi_{mn}$, $\psi_{mn}\in[0,2\pi[$, correspond to random phases. The forcing amplitude is

[TABLE]

where $N=8$ in our case. The largest forcing mode ($k_{N}=2\pi N/L$) defines the energy-carrying scale and the characteristic size of the largest magnetic eddies $l_{0}=2\pi/k_{N}$.

4.3.2 Numerical efficiency and parallel scaling

The turbulence cascade is formed as the plasma relaxes from the non-equilibrium state (Fig. 4, bottom row; see, e.g., Nättilä & Beloborodov 2021 for a more in-depth physical interpretation). We use the turbulence simulation setup to probe the numerical performance and parallel scaling of the code.

This setup is advantageous not only because it is a real physical problem, but also because we can probe both the load-balanced and load-imbalanced performance of the PIC code. Initially the simulations are in the load-balanced state (computing time being dominated by particle push time), but after the cascade forms we observe a strong numerical imbalance between different ranks because of the strong fluctuations in particle number densities (computing time dominated by communication routines). No dynamic load-balancing algorithm is used. Here we mainly focus on these transition scales of around $\sim 10^{3}$ processors where both of these regimes can be simulated, leaving extremely large simulations with $\gtrsim 10^{5}$ processors and dynamic load balancing for future work. The scaling measurement reported here is performed in the Kebnekaise computing cluster.888 Kebnekaise is a Lenovo cluster with $19288$ cores (mainly Intel Xeon E5-2690v4 Broadwell CPUs with a $2.6\leavevmode\nobreak\ \mathrm{GHz}$ base frequency) and $64\leavevmode\nobreak\ \mathrm{GB}$ of memory per each NUMA island in a node, and has FDR Infiniband switches for communication.

We use the second-order FDTD solver (FDTD2), Boris pusher (BorisPusher), linear interpolation (LinearInterpolator), and first-order ZigZag scheme (ZigZag) to perform the simulations. Calculations are done using the double floating-point precision.

Good numerical performance is obtained for the PIC code. In an ideal completely load-balanced state the average total particle update time per processor is around $0.7\leavevmode\nobreak\ \mu\mathrm{s}$ (see Fig. 5 around $n\lesssim 200$). In the strongly imbalanced simulation state (dominated mainly by inter-rank particle communications) the particle update time per processor is around $1.1\leavevmode\nobreak\ \mu\mathrm{s}$ when the turbulent cascade has developed (Fig. 5 for $n\gtrsim 4000$). The decrease in speed is caused by an increase in the evaluation time of routines communicating the fields and particles to the neighboring ranks.

We also monitor all the separate routines (see Table 1) independently to probe the numerical performance of the code in more detail. We note, however, that this monitoring is only done for one particular rank so the component-wise evaluation times reported here might fluctuate between different ranks experiencing differing loads. However, the numbers are indicative of where most of the evaluation time is spent. As seen in Fig. 5, in the ideal case most of the evaluation time (around $0.3\leavevmode\nobreak\ \mu\mathrm{s}$ per particle per processor) is spent in the interpolation routine (interp_em; pink solid line) where $\bm{\mathit{\hat{B}}}$ and $\bm{\mathit{\hat{E}}}$ are computed at the particle locations between the grid points. The large cost of the interpolation step is caused by frequent cache misses. The interpolation routine consists of many random access operations of the field arrays because particles are not sorted by their location; this prevents the processor cache prefetching from correctly predicting the incoming requested values. The second most expensive routine is the current calculation (comp_curr; green dashed line) with a typical evaluation time of $0.06\leavevmode\nobreak\ \mu\mathrm{s}$ per particle per processor. It has a similar problem with the cache memory because writing the current to the array is preformed with an unpredictable access pattern. When communication costs increase, the MPI messaging routines start to be as costly as the current deposition with an average time of about $0.1\leavevmode\nobreak\ \mu\mathrm{s}$ per particle per processor.

The code shows good weak (Fig. 6) and strong scaling (Fig. 7) when the simulations transform from cpu-bound to communication-bound. Weak scaling proceeds as $1.0:1.10:1.40:1.93:2.09$ for every quadrupling of the number of processors $16:64:256:1024:4096$. Strong scaling proceeds as $1.0:1.18:1.48:2.36:1.93$ against an ideal scaling for every doubling of the number of processors $256:512:1024:2048:4096$. The increase in the evaluation time when the number of processors is increased can indeed be attributed to the increase in the evaluation time of the communication routines (field updates: mpi_b1, mpi_b2, mpi_e1, mpi_cur, clear_vir_cur; and particle updates: mpi_prtcls, unpack_vir_prtcls, del_vir_prtcls). A similar test performed in the load-balanced state yields almost ideal scaling up to $\sim 10^{3}$ processors (scaling of roughly $1.10$ of the ideal) with a total particle update time of about $0.7\leavevmode\nobreak\ \mu\textrm{s}$.

Finally, we note that the extra costs originating from the intra-rank updates of the tile boundaries are negligible. This extra cost is introduced because of the patch-based grid infrastructure where the smallest unit of calculation is one tile ($10\times 10$ lattice with about $10^{4}$ particles per species per tile in this particular case). Therefore, some communication is always needed to keep these individual units of calculation in sync (field update routines: upd_bc0, upd_bc1, upd_bc2, cur_exchange; and particle update routines: check_outg_prtcls, pack_outg_prtcls, get_inc_prtcls, del_trnsfrd_prtcls), even in systems completely relying on shared-memory parallelism. These routines are, however, typically $10$–$100$ times cheaper than the costs of the field interpolation or MPI communication tasks.

5 Discussion

5.1 Computational advantages

runko is a modern numerical toolbox tailored for astrophysical plasma simulations. The framework is designed using recent computer science methods that rely on multiple levels of code abstraction. These in turn, help to create a general and easily extensible code framework. Furthermore, in order to encourage all kinds of use, runko is provided as an open-source software for the community.

The code is written as a hybrid C++/Python program. The low-level C++ kernels are currently designed to conform to a C++$14$ standard with the benefit of a plethora of modern computational methods. For example, we heavily rely on recent C++ template metaprogramming features to simplify the design of the code. This ensures a truly modern numerical code that will not be outdated immediately upon publication.

In addition to the underlying numerically efficient C++ implementation, all the classes have a Python interface. This allows the user to take full advantage of all the benefits of using a high-level language without sacrificing the actual run-time speed of the code. In this way, many of simulation setups can be quickly prototyped in a laptop and then scaled up to supercomputer platforms for the actual production runs. Furthermore, it allows designing and implementing very complex initial conditions for simulations because the initial routines can be created using Python.

The modular design allows runko to function as a fully agnostic framework for various different modeling formulations of simulating plasma with computers. The physics is implemented via the modules, which in practice often solve some partial differential equations or propagate Lagrangian particle species. This means that we are not locked into one predefined theoretical formalism of describing the evolution of the plasma. The high degree of modularity also allows the user to easily change the numerical solvers and algorithms depending on the problem at hand.

The implementation relies on modern C++ features such as polymorpishm and template metaprogramming to simplify the code and make it easily extensible. Different solvers inherit their standard interface from an abstract base class implementation so that introduction and implementation of new physical algorithms is straightforward. User-definable dimensionality of different solvers is handled by template specializations. Additionally, most of the solvers are agnostic to the actual floating-point accuracy (i.e., single- or double-length floating-point precision) because the working precision is typically also specified as a template variable.

Another technical advantage of the framework is the patch-based domain decomposition strategy that is built into the grid infrastructure routines. As shown by the performance analysis, the cost of the additional updates needed to keep the individual tiles in sync are negligible in comparison to the evaluation costs of the physics in modules. Having a larger number of moderately sized tiles is found to be faster than partitioning the domain with fewer but larger tiles.999For one particular test we partitioned a $3D$ PIC shock simulation grid of $7680\times 120^{2}$ with tiles of sizes $8^{3}$, $15^{3}$, or $30^{3}$ (keeping the number of processors fixed at 1280 cores) and found that the absolute wall-clock time per lap varies as $1.34$, $0.66$, and $1.22\mathrm{s}$. Smaller but more numerous tiles are found to have more expensive communication routines, whereas larger tiles have worse solver performance. The optimal size in this particular case is for tiles somewhere around $15^{3}$.

The benefits of this system are that it automatically introduces some level of load-balancing into the simulations (each rank has $\sim 10^{2}$ tiles so load imbalances at the level of one tile are effectively smoothed out), it helps in data cache locality (array sizes and particle containers remain small and easily movable in memory), and it simplifies the global communication routines (since these communication routines are anyway needed). Another benefit is the possibility to optimize the MPI domain boundaries with more complex edge shapes. An interesting example of this is the honeycomb-like tile memory configuration (see Fig. 4, top row). We also note that a similar decomposition technique currently seems to be the state-of-the-art technical choice of many other recent high-performance PIC codes, such as VPIC (Bowers et al. 2008) and Smilei (Derouillat et al. 2018).

5.2 Future directions

The presented framework offers a unique possibility to start exploring and experimenting with multiphysics simulations. These are new kinds of numerical simulations that evolve multiple physical formalisms simultaneously, or in an adaptive fashion selecting the most accurate (or relevant) solver on the fly.

These heterogeneous simulations will enable us to simulate larger, and therefore more realistic systems. In the future this can enable, for example, novel plasma simulations where the majority of the simulation domain is modeled with a force-free electrodynamic description of the plasma (Komissarov 2002). However, in some small regions with extreme conditions (e.g., current sheets, shocks) the simulation can be adaptively refined (on the fly) to the fully kinetic description that enables a realistic in situ modeling of particle acceleration. Another possibility could be to use a multiple-domain–multiple-physics approach, where, for example, some fixed part of the simulation is described with a kinetic description, whereas some other part is modeled with MHD. Here the domains could be divided based on some strong density gradient in the problem setup, like those found in systems with diluted nonthermal particles on top of a denser fluid-like plasma; physically these can model astrophysical systems like accretion disk coronas.

Another important focus is the code performance and scalability. In practice, this mostly means minimizing the global communication costs. One possibility of decreasing the communication costs is the use of more complex hybrid parallelization schemes. As an example, using the patch-based domain super-decomposition encourages the use of a task-based hybrid parallelization scheme, for example similar to the Dispatch framework (Nordlund et al. 2018). Since most of the updates of the tiles are numerically independent, these operations can be easily performed in parallel with shared-memory parallelization strategies. This, in turn, allows the number of tiles to be increased per MPI rank, which acts as an extra buffer against sudden load-imbalance fluctuations. This strategy also allows interleaved global (non-blocking) communications and normal solver updates to be performed simultaneously. This allows us to hide almost all the (MPI) communications since we can prioritize the evaluation of tiles such that the boundary tiles are always computed first and then sent to the neighbors, whereas the calculation of the inner tiles are continued independently by other threads while waiting for the halo tiles. We are also currently experimenting with dynamical load balancing methods where the ownership of the tiles changes during the simulation depending on the computational load.

Finally, since Python is becoming the new lingua franca of astronomy and astrophysics, it is important to note that problem initialization in runko can be done purely via Python scripts. This means that new users are immediately familiar with the basic approach of the code, and that initialization of very complex simulation setups can be done by a high-level language. We hope that this will accelerate the use of plasma simulations to study many yet unresolved astrophysical problems.

6 Summary

We started by reviewing the kinetic plasma theory and the numerically solved Vlasov–Maxwell system of equations (Sect. 2). The numerical implementation of the framework is discussed in Sect. 3. We focused in particular on the design of the code, and introduced the different numerical components of the framework. First, we discussed the corgi grid infrastructure and the related patch-based domain super-decomposition strategy where the computational grid is partitioned into separate tiles (Sect. 3.1). Second, we presented the physical building blocks of the code, the modules (Sect. 3.2). Different physical formulations of modeling the evolution of the plasma and electromagnetic fields are encapsulated in different modules. Third, each module can contain several numerical algorithms called solvers that are short C++ classes that operate and modify the content of different tiles (Sect. 3.3). Finally, these C++ solvers are applied to the tiles and operated by Python scripts called drivers, ensuring easy of use (Sect. 3.4).

As our first task, we implemented a new particle-in-cell module into the framework. The implementation is discussed in detail in Sects. 2.3.1 and 2.3.2. The PIC code demonstrates good numerical performance with a mean particle push time per processor of about $0.7\leavevmode\nobreak\ \mu\mathrm{s}$ in an ideal load-balanced state and $1.1\leavevmode\nobreak\ \mu\mathrm{s}$ in a strongly load-imbalanced state. Furthermore, the code is shown to scale well up to $\sim 10\,000$ processors.

We showcased the use of this PIC code by three different complex astrophysically motivated plasma simulations. First, we simulated the linear growth and nonlinear saturation of the Weibel instability in a relativistically streaming pair plasma. Second, we performed a mock simulation of a current-filamentation-mediated relativistic collisionless shock in full 3D. Last, we tested the code performance by modeling a decaying relativistic kinetic turbulence in magnetized pair plasma. These simulations demonstrated the ease of use, extendability, and flexibility of the code. All the setups are commonly used to study astrophysical phenomena, and can therefore serve as a useful starting point for new users.

Acknowledgments

JN would like to thank Andrei Beloborodov, Axel Brandenburg, Luca Comisso, Troels Haugbølle, Dhruba Mitra, Åke Nordlund, Lorenzo Sironi, Anatoly Spitkovsky, and Alexandra Veledina for stimulating discussions on numerics and plasma simulations. We would also like to thank the referee for many helpful comments and suggestions. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC and HPC2N.

Appendix A Non-dimensionalization of the Vlasov–Maxwell equations

Following Jackson (1975), the Maxwell equations in an arbitrary system of units can be written as

[TABLE]

where $k_{1}$ and $k_{2}$ are constants that define the unit system. For a Gaussian system we have $k_{1}=1$ and $k_{2}=c$, whereas for rationalized MKSA $k_{1}=1/4\pi\epsilon_{0}$ and $k_{2}=1$. Here $\epsilon_{0}$ is the vacuum permittivity (and $\mu_{0}$ is the vacuum permeability found from the relation $c^{2}=1/(\epsilon_{0}\mu_{0})$).

In this appendix we present the unit system in use (originally by Buneman et al. 1993, and described in detail also in the Tristan-MP user-manual by A. Spitkovsky). We select the Gaussian system here by setting $k_{1}=1$ and $k_{2}=c$. Our Maxwell equations are then simplified to

[TABLE]

and the fields appear symmetrical as $\bm{\mathit{E}}\leftrightarrow\bm{\mathit{B}}$. The Lorentz force in this unit system is

[TABLE]

We next normalize all the variables with some fiducial values. The most peculiar of these is our selection of distance and time scalings: we express the distance in units of the grid spacing, $\bm{\mathit{x}}=\bm{\mathit{\hat{x}}}\Delta x$, and time in units of the time step, $t=\hat{t}\Delta t$. The coordinate-velocity is then $\bm{\mathit{v}}=\bm{\mathit{\hat{v}}}\Delta x/\Delta t=\bm{\mathit{\hat{v}}}c/\hat{c}$ and the four-velocity is $\bm{\mathit{u}}=\bm{\mathit{\hat{u}}}c/\hat{c}=\bm{\mathit{\hat{v}}}\gamma(\bm{\mathit{v}})c/\hat{c}$. The fields are scaled with a fiducial field $B_{0}$ such that $\bm{\mathit{E}}=\bm{\mathit{\hat{E}}}B_{0}$ and $\bm{\mathit{B}}=\bm{\mathit{\hat{B}}}B_{0}$. Similarly, the charge, mass, and current density are also presented such that $q=\hat{q}q_{0}$, $m=\hat{m}m_{0}$, and $\bm{\mathit{J}}=J_{0}\bm{\mathit{\hat{J}}}$, respectively. This way, all of our numerical quantities are denoted with a hat.

In numerical (code) units the equations that are being solved are then

[TABLE]

where in the last steps we define

[TABLE]

so that Eqs. (33)–(35) appear unitless.

From Eqs. (36) and (37) we obtain a connection between the grid spacing, particle mass, and charge. We note first that the current density is

[TABLE]

i.e., $J_{0}=q_{0}/\Delta x^{2}\Delta t$. This implies that

[TABLE]

i.e.,

[TABLE]

Assuming that one computational macroparticle contains $N$ electrons (or positrons) we can write

[TABLE]

where we have set $m_{0}=m_{\mathrm{e}}$ ($m_{\mathrm{e}}$ is the electron rest mass) and $q_{0}=e$ ($e$ is the electron charge), and where $r=e^{2}/m_{\mathrm{e}}c^{2}\approx 2.82\times 10^{-13}\,\mathrm{cm}$ is the classical electron radius. This enables us to express the field scaling, Eq. (36), in physical (Gaussian) units

[TABLE]

The conversions from code units to Gaussian units can then be performed by selecting a length scale, $\Delta x$, as

[TABLE]

We next discuss some derived plasma quantities in the simulation. The total relativistic plasma frequency is given as

[TABLE]

where $\omega_{\mathrm{p},-}$ and $\omega_{\mathrm{p},+}$ are the plasma frequencies of electrons and ions (or electrons and positrons), respectively, and where we have assumed $n=n_{-}=n_{+}$. For electron-ion plasma $m_{-}=1$ and $m_{+}=1836$ (or some smaller reduced proton mass), whereas for electron-positron pair plasma $m_{-}=1$, $m_{+}=1$, and $1+m_{-}/m_{+}=2$. Time dilatation from a relativistic bulk-motion (if any) is corrected by the mean Lorentz factor of that flow $\langle\gamma\rangle$. The counter-part in the code units is

[TABLE]

where $N_{ppc}$ is the number of particles per cell per species (again assuming charge equilibrium, $N_{ppc} =n_{-}=n_{+}$) and in the last step we take into account that $|\hat{q}|/\hat{m}=1$. Initial numerical plasma oscillation frequency is obtained by specifying the skin-depth resolution, $\lambda_{\mathrm{p}}=c/\omega_{\mathrm{p}}=\mathcal{R}_{\mathrm{p}}\Delta x$, and remembering that $c=\hat{c}\Delta x/\Delta t$ and $\omega_{\mathrm{p}}=\hat{\omega_{\mathrm{p}}}/\Delta t$, so that we obtain

[TABLE]

We can then fix the value of charge $|\hat{q}|$ (and $\hat{m}$) such that the initial plasma frequency for the electrons in the beginning of the simulation is (i.e., requiring $\Delta t=\hat{\omega}_{\mathrm{p},0}^{-1})$

[TABLE]

Similarly, the relativistic species-specific cyclotron frequency is

[TABLE]

which in code units is simply

[TABLE]

where $\hat{m}\frac{m_{+}}{m_{-}}$ is the numerical mass of the particle (i.e., $1\hat{m}$ for pair plasma or $1836\hat{m}$ for protons). We note that the critical electron cyclotron frequency for $B=B_{0}$ (i.e., $\bm{\mathit{\hat{B}}}=1$) corresponds to $\omega_{B,0}\Delta t=1/\hat{c}$. Both Eqs. (48) and (52) can be translated to physical units by simply inserting numerical values to the time step $\Delta t^{-1} =c/\hat{c}\Delta x$. By comparing Eqs. (47) vs. (48) and Eqs. (51) vs. (52) we also see that numerical quantities are conveniently obtained by just replacing the physical quantities with the code variables.

Finally, the total plasma magnetization is given as

[TABLE]

This demonstrates the usefulness of the unit system. Numerically we perform fewer floating-point operations because formulas appear unitless, our typical values are of the same order of magnitude reducing floating-point round-off errors, and physical quantities are easily transformed to the code unit system by replacing variables with code quantities (e.g., $x\rightarrow\hat{x}$ and $n=N_{ppc}$) and dropping $4\pi$ factors from the Maxwell equations.

Appendix B Time and space staggering with finite-difference method

Numerical charge conservation and divergence-free $\bm{\mathit{B}}$ field preservation follow from the use of a charge conserving current deposition strategy and the staggered Yee lattice (Yee 1966). See Esirkepov (2001) for more in-depth discussion of the topic.

Let us define discrete versions of the $\nabla$ operator as

[TABLE]

and

[TABLE]

where $\Delta x$, $\Delta y$, and $\Delta z$ are the grid spacings. These operators have the following properties

[TABLE]

where $\hat{\Delta}$ is the discrete Poisson operator in central differences such that

[TABLE]

Let us stagger our variables in space and time such that

[TABLE]

Maxwell’s equations for $\bm{\mathit{E}}$ and $\bm{\mathit{B}}$ are given as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The discrete form of the charge conservation then follows by operating with $\nabla^{+}\cdot$ on Eq. (64) and substituting (66) into the equation that follows. We then obtain

[TABLE]

Similarly, operating on Eq. (65) with $\nabla^{-}$ we obtain the divergence-free nature of $\bm{\mathit{B}}$ field

[TABLE]

Thus, if

[TABLE]

and

[TABLE]

at $t=0$, then the divergence of $\bm{\mathit{E}}$ is always equal to the charge density and $\bm{\mathit{B}}$ will retain its divergence-free nature up to machine precision during the simulation.

Appendix C Relativistic Boris scheme

We present the full relativistic Boris scheme here that updates particle four-velocity one time step forward. This can be found as Pusher solver named BorisPusher in the pic module. In the beginning we have a particle with a four-velocity $\bm{\mathit{u}}^{n}$ in units of $c$. The Lorentz update is done via velocities in units of $\hat{c}$ since this is the effective speed of light experienced by the electromagnetic fields. To differentiate between these two units we use primed quantities to mark velocities in units of $\hat{c}$.

The scheme begins with a first half of electric field acceleration as

[TABLE]

where

[TABLE]

and $m_{\pm}=m/m_{\mathrm{e}}$ is the mass of the particle in units of electron rest mass $m_{\mathrm{e}}$. This four-velocity in units of $\hat{c}$ corresponds to a Lorentz factor of

[TABLE]

We then proceed with a first half of magnetic rotation

[TABLE]

where

[TABLE]

As a final step we combine the second half of magnetic rotation and second half of the electric field acceleration as

[TABLE]

Final particle four-velocity in units of $c$ is then

[TABLE]

summarizing a recipe for updating the particle velocity from an initial time $\bm{\mathit{u}}^{n}$ to $\bm{\mathit{u}}^{n+1}$ (i.e., one time step of $\Delta t$ forward).

Appendix D Current deposition

We summarize here the charge-conserving ZigZag current deposition scheme (Umeda et al. 2003) implemented as the ZigZag solver in runko. First, we compute the relay point $\bm{\mathit{x}}_{\mathrm{r}}$ between previous particle position $\bm{\mathit{x}}_{1}$ (and closest grid index $\bm{\mathit{i}}_{1}=(i_{1},j_{1},k_{1})$) and the current position $\bm{\mathit{x}}_{2}$ (and closest grid index $\bm{\mathit{i}}_{2}$) as

[TABLE]

Here we assume element-wise $\mathrm{min}$ and $\mathrm{max}$ operators as $\mathrm{min}\{\bm{\mathit{a}},\bm{\mathit{b}}\}\equiv(\mathrm{min}\{a_{x},b_{x}\},\mathrm{min}\{a_{y},b_{y}\},\mathrm{min}\{a_{z},b_{z}\})$.

The current fluxes are then obtained as

[TABLE]

We assume first-order shape functions for the particles. This is related to particle weights defined at the mid points $(\bm{\mathit{x}}_{\mathrm{r}}+\bm{\mathit{x}}_{1})/2$ and $(\bm{\mathit{x}}_{2}+\bm{\mathit{x}}_{r})/2$ as

[TABLE]

Depending on the dimensionality, the weights are described as $\bm{\mathit{W}}_{1,2}=(W_{1,2;x})$ for $D=1$, $\bm{\mathit{W}}_{1,2}=(W_{1,2;x},W_{1,2;y})$ for $D=2$, or $\bm{\mathit{W}}_{1,2}=(W_{1,2;x},W_{1,2;y},W_{1,2;z})$ for $D=3$. Finally, the current density $J=FW$ needs to be injected into adjacent grid points following a charge-conserving deposition strategy. This can be implemented following Umeda et al. (2003).