# ESVM: an open-source finite volume Electrostatic Vlasov-Maxwell code

Michaël J. Touati¹, ², ³

¹ Department of Electrical Engineering, University of California Los Angeles, Los Angeles, CA 90095, USA ² Group of Lasers and Plasmas, IPFN, IST, Universidade de Lisboa, Lisbon, Portugal ³ Centro de Láseres Pulsados de Salamanca (CLPU), Edificio M5, Parque Científico, C/ Adaja 8, 37185 Villamayor, Salamanca, Spain (current affiliation)

DOI: 10.21105/joss.03618

## Software

- Review ⚗

- Repository ⚗

- Archive ⚗

Editor: Daniel S. Katz ⚗

## Reviewers:

- @tclune

- @TomGoffrey

- @rouson

Submitted: 11 August 2021

Published: 29 November 2021

## License

Authors of papers retain copyright and release the work under a Creative Commons Attribution 4.0 International License (CC BY 4.0).

## Summary

A plasma is a set of charged particles consisting of electrons and ionized atoms whose quantity is sufficiently large to behave collectively through the long-distance electromagnetic fields they produce. It is thought that more than 99.9% of visible matter in the universe is in the plasma state. In a collisionless plasma consisting in an ionized gas composed of electrons moving in-between much heavier ions, any electrostatic field is rapidly screened by the plasma electrons over the Debye screening distance (Debye & Hückel, 1923). When the number of electrons in these Debye spheres can be assumed to be infinite, the plasma electron population is correctly described by the Vlasov equation (Vlasov, 1938) that neglects all correlations between particles such as the binary Coulomb collisions between them. In addition to being simple, the resulting Vlasov-Maxwell set of equations is extremely rich in physics and has many applications ranging from astrophysics and theoretical plasma physics to intense laser-matter interaction experiments. ESVM (ElectroStatic Vlasov-Maxwell) is a Vlasov-Maxwell Fortran 95 standard-compliant code, parallelized with OpenMP and using Python 3 for post-processing, that allows for the study of these collisionless plasmas. Many finite volume advection schemes (Godunov, 1959) are implemented in order to discretize the Vlasov equation, namely:

- the donor-cell scheme, i.e., the downwind / upwind scheme (Courant et al., 1952) depending on the advection direction in each phase-space cell,

- the Lax-Wendroff scheme (Lax & Wendroff, 1960),

- the Fromm scheme (Fromm, 1968),

- the Beam-Warming scheme (Beam & Warming, 1976),

- the Van Leer scheme (Van Leer, 1977),

- the minmod scheme (Roe, 1986),

- the superbee scheme (Roe, 1986), and

- two Monotonic Upwind-centered Scheme for Conservation Laws (MUSCL) (van Leer, 1979) schemes MUSCL2 (Crouseilles & Filbet, 2004) and MUSCL1 (Duclous et al., 2009).

Unlike the linear second order Lax-Wendroff, Fromm, and Beam-Warming schemes, the non-linear second order minmod, superbee, Van Leer, and MUSCL schemes make use of a Total Variation Diminishing (TVD) non-linear flux limiter with the price of becoming a first order scheme in some phase-space cells to limit the numerical oscillations. The donor-cell scheme is a first order method and has the pros of limiting such eventual oscillations but the cons of being numerically less consistent and more diffusive. In ESVM, the discretized Vlasov equation is coupled with the self-consistent Maxwell-Gauss equation or equivalently with the Maxwell-Ampere equation with Maxwell-Gauss equation computed at the first time step, only.
While the second order Maxwell-Gauss solver needs a computationally expensive inversion of a tridiagonal matrix for the computation of the Poisson equation, the Maxwell-Ampere equation solver makes use of a faster first order numerical scheme (in time). Both absorbing and periodic boundary conditions for both the particles and the fields are implemented. Python scripts, using the Matplotlib and Numpy packages, are provided to automatically extract and plot the stored simulation results. The simulation parameters are described in the *input deck* and they can be modified without having to recompile the code. Compilation rules can be modified in the *makefile* depending on the user's compiler preferences. Classical plasma physics academic case simulations that need less than one CPU × hour each, tools for testing the compilation of the code, and tools for checking the simulation parameters are provided.

## Statement of need

ESVM has been developed in order to adapt simulations to specific plasma physics problems by choosing the more adequate finite volume numerical advection scheme in order to compute the Vlasov equation phase-space advection derivatives and to choose between computing the Maxwell-Gauss equation or the Maxwell-Ampere equation with Maxwell-Gauss equation computed at the first time step, only. The code aims at beeing used by the open-source highly parallel computing (HPC) plasma physics community, ranging from under or post-graduate students to teachers and researchers who usually use Particle-In-Cell (PIC) codes (Dawson, 1962) to study collisionless plasmas. Indeed, the PIC method may prohibit the study of plasma physical processes on large time scales and/or for very dense collisionless plasmas due to the statistical and numerical fluctuations of the computed quantities imposed by the use of a finite number of particles. Also, plasma instabilities naturally develop in PIC codes, seeded by the available fluctuations spatial spectrum k-vector for which the instability growth rate is maximum and some small amplitude plasma physical processes may be hidden under the fluctuations level. Compared to the many open source PIC codes such as Smilei (Derouillat et al., 2018) and semi-Lagrangian codes such as vmf90 (de Buyl, 2014), there are no open source finite volume Vlasov-Maxwell codes in the literature that are not based on an expansion method such as OSHUN (Tzoufras et al., 2011), AMORE (Touati et al., 2014), or Vlapy (Joglekar & Levy, 2020). Finally, since the Vlasov equation is a conservation equation of the probable number of particles in the phase-space, using a finite volume method in order to compute the Vlasov equation has the advantage of allowing for the use of numerical schemes that are numerically flux conserving and/or that ensure the distribution function positivity compared to other numerical methods. ESVM has already been used during courses for under and post-graduate students about the "numerical tools for laser-plasma interaction Physics" and it is currently used for theoretical Plasma Physics investigations.

## Future work

The author plans, in the near future, to:

1. provide another plasma physics academic simulation about one non-linear BGK electron plasma wave, from the name of its finder I. B. Bernstein, J. M. Greene, and M. D. Kruskal (Bernstein et al., 1957)

2. provide another plasma physics academic simulation about the echo of two plasma electron waves (Gould et al., 1967)

3. implement non-eqally spaced phase-space cells

4. implement high order Weighted Essentially Non-Oscillatory (WENO) advection schemes (Liu et al., 1994)

5. compute the plasma ion Vlasov equation to allow for the ions to be mobile
6. store the simulation results in hdf5 files instead of text files

7. implement MPI parallelization

8. implement vectorization

9. extend the code to the relativistic regime: ESVM ⇒ RESVM for open source Relativistic ElectroStatic Vlasov-Maxwell code

10. implement a BGK collision operator, from the name of its finder P. L. Bhatnagar, E. P. Gross, and M. Krook (Bhatnagar et al., 1954)

11. extend the code to 1D-2V and 1D-3V phase-space electrostatic plasma simulations

12. implement the Landau (Landau, 1937) and Belaiev-Budker (Belaiev & Budker, 1956) relativistic collision operators using the Rosenbluth potentials (Rosenbluth et al., 1957) and their relativistic Braams-Karney extension (Braams & Karney, 1987): (R)ESVM ⇒ (R)ESVFPM for open source (Relativistic) ElectroStatic Vlasov-Fokker-Planck-Maxwell code

13. extend the code to electrostatic 2D-1V, 2D-2V, and 2D-3V phase-space plasma simulations: (R)ESV(FP)M ⇒ (R)ESV(FP)M2 for open source (Relativistic) ElectroStatic Vlasov-(Fokker-Planck)-Maxwell in 2D

14. extend the code with the second order finite difference Yee scheme (Yee, 1966) to electromagnetic 2D-1V, 2D-2V, and 2D-3V phase-space plasma simulations: (R)ESV(FP)M(2) ⇒ (R)EMV(FP)M(2) for open source (Relativistic) ElectroMagnetic Vlasov-(Fokker-Planck)-Maxwell (in 2D)

15. implement the Perfectly Matched Layer (PML) technique (Berenger, 1994) to absorb the electromagnetic fields at the spatial simulation box boundaries

16. deploy the code to GPU architectures.

## References

Beam, R. M., & Warming, R. F. (1976). An implicit finite-difference algorithm for hyperbolic systems in conservation-law form. *Journal of Computational Physics*, **22**(1), 87–110.
https://doi.org/10.1016/0021-9991(76)90110-8

Belaiev, S. T., & Budker, G. I. (1956). *Doklady Akademii Nauk SSSR*, **107**.

Berenger, J.-P. (1994). A perfectly matched layer for the absorption of electromagnetic waves. *Journal of Computational Physics*, **114**(2), 185–200. https://doi.org/10.1006/jcph.1994.1159

Bernstein, I. B., Greene, J. M., & Kruskal, M. D. (1957). Exact nonlinear plasma oscillations. *Physical Review*, **108**, 546–550. https://doi.org/10.1103/PhysRev.108.546

Bhatnagar, P. L., Gross, E. P., & Krook, M. (1954). A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. *Physical Review*, **94**, 511–525. https://doi.org/10.1103/PhysRev.94.511

Braams, B. J., & Karney, C. F. F. (1987). Differential form of the collision integral for a relativistic plasma. *Physical Review Letter*, **59**, 1817–1820. https://doi.org/10.1103/PhysRevLett.59.1817

Courant, R., Isaacson, E., & Rees, M. (1952). On the solution of nonlinear hyperbolic differential equations by finite differences. *Communications on Pure and Applied Mathematics*, **5**(3), 243–255. https://doi.org/10.1002/cpa.3160050303

Crouseilles, N., & Filbet, F. (2004). Numerical approximation of collisional plasmas by high order methods. *Journal of Computational Physics*, **201**(2), 546–572. https://doi.org/10.1016/j.jcp.2004.06.007

Dawson, J. (1962). One-dimensional plasma model. *The Physics of Fluids*, **5**(4), 445–459.
https://doi.org/10.1063/1.1706638
de Buyl, P. (2014). The vmf90 program for the numerical resolution of the vlasov equation for mean-field systems. *Computer Physics Communications*, **185**(6), 1822-1827. https://doi.org/10.1016/j.cpc.2014.03.004

Debye, P., & Hückel, E. (1923). Zur theorie der elektrolyte. I. Gefrierpunktserniedrigung und verwandte erscheinungen. *Zeitschrift für Physik*, **24**, 305–324.

Derouillat, J., Beck, A., Pérez, F., Vinci, T., Chiaramello, M., Grassi, A., Flé, M., Bouchard, G., Plotnikov, I., Aunai, N., Dargent, J., Riconda, C., & Grech, M. (2018). Smilei : A collaborative, open-source, multi-purpose particle-in-cell code for plasma simulation. *Computer Physics Communications*, **222**, 351–373. https://doi.org/10.1016/j.cpc.2017.09.024

Duclous, R., Dubroca, B., Filbet, F., & Tikhonchuk, V. (2009). High order resolution of the maxwell-fokker-planck-landau model intended for ICF applications. *Journal of Computational Physics*, **228**(14), 5072–5100. https://doi.org/10.1016/j.jcp.2009.04.005

Fromm, J. E. (1968). A method for reducing dispersion in convective difference schemes. *Journal of Computational Physics*, **3**(2), 176–189. https://doi.org/10.1016/0021-9991(68)90015-6

Godunov, S. K. (1959). Eine differenzenmethode für die näherungsberechnung unstetiger lösungen der hydrodynamischen gleichungen. *Matematicheskii Sbornik. Novaya Seriya*, **47**, 271–306.

Gould, R. W., O'Neil, T. M., & Malmberg, J. H. (1967). Plasma wave echo. *Physical Review Letter*, **19**, 219-222. https://doi.org/10.1103/PhysRevLett.19.219

Joglekar, A. S., & Levy, M. C. (2020). VlaPy: A python package for eulerian vlasov-poisson-fokker-planck simulations. *Journal of Open Source Software*, **5**(53), 2182. https://doi.org/10.21105/joss.02182

Landau, L. D. (1937). *Journal of Experimental and Theoretical Physics*, **7**, 203.

Lax, P., & Wendroff, B. (1960). Systems of conservation laws. *Communications on Pure and Applied Mathematics*, **13**(2), 217–237. https://doi.org/10.1002/cpa.3160130205

Liu, X.-D., Osher, S., & Chan, T. (1994). Weighted essentially non-oscillatory schemes. *Journal of Computational Physics*, **115**(1), 200–212. https://doi.org/10.1006/jcph.1994.1187

Roe, P. L. (1986). Characteristic-based schemes for the euler equations. *Annual Review of Fluid Mechanics*, **18**(1), 337–365. https://doi.org/10.1146/annurev.fl.18.010186.002005

Rosenbluth, M. N., MacDonald, W. M., & Judd, D. L. (1957). Fokker-planck equation for an inverse-square force. *Physical Review*, **107**, 1–6. https://doi.org/10.1103/PhysRev.107.1

Touati, M., Feugeas, J.-L., Nicolaï, P., Santos, J. J., Gremillet, L., & Tikhonchuk, V. T. (2014). A reduced model for relativistic electron beam transport in solids and dense plasmas. *New Journal of Physics*, **16**(7), 073014. https://doi.org/10.1088/1367-2630/16/7/073014

Tzoufras, M., Bell, A. R., Norreys, P. A., & Tsung, F. S. (2011). A vlasov-fokker-planck code for high energy density physics. *Journal of Computational Physics*, **230**(17), 6475–6494. https://doi.org/10.1016/j.jcp.2011.04.034

van Leer, B. (1979). Towards the ultimate conservative difference scheme. V. A second-order sequel to godunov's method. *Journal of Computational Physics*, **32**(1), 101–136. https://doi.org/10.1016/0021-9991(79)90145-1

Van Leer, B. (1977). Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow. *Journal of Computational Physics*, **23**(3), 263–275. https://doi.org/10.1016/0021-9991(77)90094-8
Vlasov, A. A. ;. (1938). *Journal of Experimental and Theoretical Physics*, **8**(3), 291.

Yee, K. (1966). *IEEE Transactions on Antennas and Propagation*, *14*(3), 302–307. https://doi.org/10.1109/TAP.1966.1138693