A massively parallel semi-Lagrangian solver for the six-dimensional Vlasov-Poisson equation

This paper presents an optimized and scalable semi-Lagrangian solver for the Vlasov-Poisson system in six-dimensional phase space. Grid-based solvers of the Vlasov equation are known to give accurate results. At the same time, these solvers are challenged by the curse of dimensionality resulting in very high memory requirements, and moreover, requiring highly efficient parallelization schemes. In this paper, we consider the 6d Vlasov-Poisson problem discretized by a split-step semi-Lagrangian scheme, using successive 1d interpolations on 1d stripes of the 6d domain. Two parallelization paradigms are compared, a remapping scheme and a classical domain decomposition approach applied to the full 6d problem. From numerical experiments, the latter approach is found to be superior in the massively parallel case in various respects. We address the challenge of artificial time step restrictions due to the decomposition of the domain by introducing a blocked one-sided communication scheme for the purely electrostatic case and a rotating mesh for the case with a constant magnetic field. In addition, we propose a pipelining scheme that enables to hide the costs for the halo communication between neighbor processes efficiently behind useful computation. Parallel scalability on up to 65k processes is demonstrated for benchmark problems on a supercomputer

A mass-conserving sparse grid combination technique with biorthogonal hierarchical basis functions for kinetic simulations

The exact numerical simulation of plasma turbulence is one of the assets and challenges in fusion research. For grid-based solvers, sufficiently fine resolutions are often unattainable due to the curse of dimensionality. The sparse grid combination technique provides the means to alleviate the curse of dimensionality for kinetic simulations. However, the hierarchical representation for the combination step with the state-of-the-art hat functions suffers from poor conservation properties and numerical instability. The present work introduces two new variants of hierarchical multiscale basis functions for use with the combination technique: the biorthogonal and full weighting bases. The new basis functions conserve the total mass and are shown to significantly increase accuracy for a finite-volume solution of constant advection. Further numerical experiments based on the combination technique applied to a semi-Lagrangian Vlasov--Poisson solver show a stabilizing effect of the new bases on the simulations

Afterlive: A performant code for Vlasov-Hybrid simulations

A parallelized implementation of the Vlasov-Hybrid method [Nunn, 1993] is presented. This method is a hybrid between a gridded Eulerian description and Lagrangian meta-particles. Unlike the Particle-in-Cell method [Dawson, 1983] which simply adds up the contribution of meta-particles, this method does a reconstruction of the distribution function $f$ in every time step for each species. This interpolation method combines meta-particles with different weights in such a way that particles with large weight do not drown out particles that represent small contributions to the phase space density. These core properties allow the use of a much larger range of macro factors and can thus represent a much larger dynamic range in phase space density. The reconstructed phase space density $f$ is used to calculate momenta of the distribution function such as the charge density $\rho$. The charge density $\rho$ is also used as input into a spectral solver that calculates the self-consistent electrostatic field which is used to update the particles for the next time-step. Afterlive (A Fourier-based Tool in the Electrostatic limit for the Rapid Low-noise Integration of the Vlasov Equation) is fully parallelized using MPI and writes output using parallel HDF5. The input to the simulation is read from a JSON description that sets the initial particle distributions as well as domain size and discretization constraints. The implementation presented here is intentionally limited to one spatial dimension and resolves one or three dimensions in velocity space. Additional spatial dimensions can be added in a straight forward way, but make runs computationally even more costly.Comment: Accepted for publication in Computer Physics Communication

Solving the Vlasov–Maxwell equations using Hamiltonian splitting

In this paper, the numerical discretizations based on Hamiltonian splitting for solving the Vlasov–Maxwell system are constructed. We reformulate the Vlasov–Maxwell system in Morrison–Marsden–Weinstein Poisson bracket form. Then the Hamiltonian of this system is split into five parts, with which five corresponding Hamiltonian subsystems are obtained. The splitting method in time is derived by composing the solutions to these five subsystems. Combining the splitting method in time with the Fourier spectral method and finite volume method in space gives the full numerical discretizations which possess good conservation for the conserved quantities including energy, momentum, charge, etc. In numerical experiments, we simulate the Landau damping, Weibel instability and Bernstein wave to verify the numerical algorithms

On the error propagation of semi-Lagrange and Fourier methods for advection problems

In this paper we study the error propagation of numerical schemes for the advection equation in the case where high precision is desired. The numerical methods considered are based on the fast Fourier transform, polynomial interpolation (semi-Lagrangian methods using a Lagrange or spline interpolation), and a discontinuous Galerkin semi-Lagrangian approach (which is conservative and has to store more than a single value per cell). We demonstrate, by carrying out numerical experiments, that the worst case error estimates given in the literature provide a good explanation for the error propagation of the interpolation-based semi-Lagrangian methods. For the discontinuous Galerkin semi-Lagrangian method, however, we find that the characteristic property of semi-Lagrangian error estimates (namely the fact that the error increases proportionally to the number of time steps) is not observed. We provide an explanation for this behavior and conduct numerical simulations that corroborate the different qualitative features of the error in the two respective types of semi-Lagrangian methods. The method based on the fast Fourier transform is exact but, due to round-off errors, susceptible to a linear increase of the error in the number of time steps. We show how to modify the Cooley--Tukey algorithm in order to obtain an error growth that is proportional to the square root of the number of time steps. Finally, we show, for a simple model, that our conclusions hold true if the advection solver is used as part of a splitting scheme.Comment: submitted to Computers & Mathematics with Application

Practical 3-splitting beyond Strang

Operator splitting is a popular divide-and-conquer strategy for solving differential equations. Typically, the right-hand side of the differential equation is split into a number of parts that can then be integrated separately. Many methods are known that split the right-hand side into two parts. This approach is limiting, however, and there are situations when 3-splitting is more natural and ultimately more advantageous. The second-order Strang operator-splitting method readily generalizes to a right-hand side splitting into any number of operators. It is arguably the most popular method for 3-splitting because of its efficiency, ease of implementation, and intuitive nature. Other 3-splitting methods exist, but they are less well-known, and evaluations of their performance in practice are scarce. We demonstrate the effectiveness of some alternative 3-split, second-order methods to Strang splitting on two problems: the reaction-diffusion Brusselator, which can be split into three parts that each have closed-form solutions, and the kinetic Vlasov--Poisson equations that is used in semi-Lagrangian plasma simulations. We find alternative second-order 3-operator-splitting methods that realize efficiency gains of 10\%--20\% over traditional Strang splitting

hyper.deal: An efficient, matrix-free finite-element library for high-dimensional partial differential equations

This work presents the efficient, matrix-free finite-element library hyper.deal for solving partial differential equations in two to six dimensions with high-order discontinuous Galerkin methods. It builds upon the low-dimensional finite-element library deal.II to create complex low-dimensional meshes and to operate on them individually. These meshes are combined via a tensor product on the fly and the library provides new special-purpose highly optimized matrix-free functions exploiting domain decomposition as well as shared memory via MPI-3.0 features. Both node-level performance analyses and strong/weak-scaling studies on up to 147,456 CPU cores confirm the efficiency of the implementation. Results of the library hyper.deal are reported for high-dimensional advection problems and for the solution of the Vlasov--Poisson equation in up to 6D phase space.Comment: 33 pages, 18 figure