VlaPy: A Python package for Eulerian Vlasov-Poisson-Fokker-Planck Simulations

The Vlasov-Poisson-Fokker-Planck system of equations is commonly used to study plasma and fluid physics in a broad set of topical environments, ranging from space physics, to laboratorycreated plasmas for fusion applications (Betti & Hurricane, 2016; Chen, Klein, & Howes, 2019; Fasoli et al., 2016; Ongena, Koch, Wolf, & Zohm, 2016). More specifically, the Vlasov-Poisson system of equations is typically employed to model collisionless plasmas. The Fokker-Planck operator can be introduced into this system to represent the effect of collisions. The primary advantage of this scheme is that instead of relying on numerical diffusion to smooth smallscale structures that arise when modeling collisionless plasmas, the Fokker-Planck operator enables a physics-based smoothing mechanism.


Statement of Need
The 1D-1V VPFP equation set solved here has been applied in research on laser-plasma interactions in the context of inertial fusion (Banks, Brunner, Berger, & Tran, 2016;Fahlen, Winjum, Grismayer, & Mori, 2009), plasma-based accelerators (Thomas, 2016), space physics (Chen et al., 2019), and fundamental plasma physics (Heninger & Morrison, 2018;Pezzi, Valentini, & Veltri, 2016). While there are VPFP software libraries which are available in academic settings, research laboratories, and industry (Banks, Brunner, Berger, Arrighi, & Tran, 2017;Joglekar et al., 2018), the community has yet to benefit from a simple-to-read, open-source Python implementation. This lack of capability is currently echoed in conversations within the PlasmaPy (PlasmaPy Community et al., 2018) community (PlasmaPy is a collection of open-source plasma physics resources). Our aim with VlaPy is to take a step towards filling this need for a research and educational tool in the open-source community.
VlaPy is intended to help students learn fundamental concepts and help researchers discover novel physics and applications in plasma physics, fluid physics, computational physics, and numerical methods. It is also designed to provide a science-accessible introduction to industry and software engineering best-practices, including unit and integrated testing, and extensible and maintainable code.
The details of the VlaPy implementation are provided in the following sections.

Equations
The Vlasov-Poisson-Fokker-Planck system can be decomposed into four components. These components, represented using normalized units, th where v th is the thermal velocity, ω p is the electron plasma frequency, m e is the electron mass, λ D is the Debye length, and e is the elementary charge. The Fourier transform operator is represented by F and the subscript to the operator indicates the dimension of the transform. In what follows, we have omitted the tilde for brevity.

Vlasov Equation
The normalized, non-relativistic (γ = 1) Vlasov equation for electrons is given by is the electron velocity distribution function.
We use operator splitting to advance the time-step (Strang, 1968). Each one of those operators is then integrated pseudo-spectrally using the following methodology.
We use the Fourier expansions of the distribution function, which are given by We first discretize f (x, v, t) = f n (x l , v j ), and then perform a Fourier expansion inx for each grid value of v.

This gives
which is substituted into the Fourier transform of the advection operator inx, as given by This process enables the decoupling ofx andv grids from the time dimension and allows us to write an Ordinary Differential Equation in time for the discretized distribution function f x n (k x , v j ). This is given by Next, we solve for the change in the plasma distribution function, integrate in time, and evaluate the integral atf x n andf x n+1 which giveŝ The E∂f /∂v term is evolved similarly usinĝ We have implemented a simple Leapfrog scheme as well as a 4th order integrator called the Position-Extended-Forest-Ruth-Like Algorithm (PEFRL) (Omelyan, Mryglod, & Folk, 2002) Tests The implementation of this equation is tested in the integrated tests section.

Poisson Equation
The normalized Poisson equation is simply Because the ion species are effectively static and form a charge-neutralizing background to the electron dynamics, we can express the Poisson equation as This is justifed by the assumption that the relevant time-scales are short compared to the time-scale associated to ion motion.
In one spatial dimension, this can be expressed as and the discretized version that is solved is

Integrated Code Testing
Unit tests are provided for this operator to validate its performance and operation under the above assumptions. These are simply unit tests against analytical solutions of integrals of periodic functions. They can be found in tests/test_fieldsolver.py.
Below, we provide an example illustration of this validation. The code is provided in notebo oks/test_poisson.ipynb Testing the Field Solver
The first of these implementations (LB) has the governing equation given by ( δf δt is the thermal velocity of the distribution. The second of these implementations (DG) has a governing equation given by is the mean velocity of the distribution and is the thermal velocity of the shifted distribution.
The second implementation is an extension of the first, and extends momentum conservation for distributions that have a non-zero mean velocity.
We discretize this backward-in-time, centered-in-space. This procedure results in the time-step scheme given by This forms a tridiagonal system of equations that can be directly inverted.

Integrated Code Testing
Unit tests are provided for this operator. They can be found in tests/test_lb.py and tests/test_dg.py. The unit tests ensure that 1. The operator does not impact a Maxwell-Boltzmann distribution already satisfying v th = v 0 .
2. The LB operator conserves number density, momentum, and energy when initialized with a zero mean velocity.
3. The DG operator conserves number density, momentum, and energy when initialized with a non-zero mean velocity.
The notebooks/test_fokker_planck.ipynb notebook contains illustrations and examples for these tests. Below, we show results from some of the tests for illustrative purposes. Testing for density, momentum, and energy conservation for non-zero mean velocity distribution with the Dougherty operator Density conservation error = 7.04e-08 Momentum conservation error = 4.23e-07 Energy conservation error = 2.55e-06 We see from the above figures that the distribution relaxes to a Maxwellian. Depending on the implementation, certain characteristics of momentum conservation are enforced or avoided.

Integrated Code Tests against Plasma Physics: Electron Plasma Waves and Landau Damping
Landau Damping is one of the most fundamental phenomena in plasma physics. An extensive review is provided in (Ryutov, 1999).
Plasmas can support electrostatic oscillations. The oscillation frequency is given by the elec-trostatic electron plasma wave (EPW) dispersion relation. When a wave of sufficiently small amplitude is driven at the resonant wave-number and frequency pairing, there is a resonant exchange of energy between the plasma and the electric field, and the electrons can damp the electric field. The damping rates, as well as the resonant frequencies, are given in (Canosa, 1973).
In the VlaPy simulation code, we have verified that the known damping rates for Landau Damping are reproduced, for a few different wave-numbers. This is shown in notebooks/la ndau_damping.ipynb.
We include validation against this phenomenon as an automated integrated test. The tests can be found in tests/test_landau_damping.py Below, we also illustrate a manual validation of this phenomenon through the fully integrated workflow of running a simulation on a local machine and sending the results to the MLFlowdriven logging mechanism. After running a properly initialized simulation, we show that the damping rate of an electron plasma wave with k = 0.3 is reproduced accurately through the UI. This can also be computed manually (please see the testing code for details).