Search

Semi-Lagrangian solvers for the Vlasov system offer noiseless solutions compared to Lagrangian particle methods and can handle larger time steps compared to Eulerian methods. In order to reduce the computational complexity of the interpolation steps, it is common to use a directional splitting. However, this typically yields the wrong angular velocity. In this paper, we analyze a semi-Lagrangian method that treats the $v \times B$ term with a rotational grid and combines this with a directional splitting for the remaining terms. We analyze the convergence properties of the scheme both analytically and numerically. The favorable numerical properties of the rotating grid solution are demonstrated for the case of ion Bernstein waves.

With the increase in computational capabilities over the last years it becomes possible to simulate more and more complex and accurate physical models. Gyrokinetic theory has been introduced in the 1960s and 1970s in the need of describing a plasma with more accurate models than fluid equations, but eliminating the complexity of the fast gyration about the magnetic field lines. Although results from current gyrokinetic computer simulations are in fair agreement with experimental results in core physics, crucial assumptions made in the derivation make it unreliable in regimes of higher fluctuations and stronger gradient, such as the tokamak edge. With our novel optimized and scalable semi-Lagrangian solver we are able to simulate ion-temperature gradient modes with the 6D kinetic model including the turbulent saturation. After thoroughly testing our simulation code against analytical computations and gyrokinetic simulations (with the gyrokinetic code GYRO), it has been possible to show first plasma properties that go beyond standard gyrokinetic simulations. This includes the explicit description of the complete perpendicular energy fluxes and the excitation of high frequency waves (around the Larmor frequency) in the nonlinear saturation phase.

In this paper, we describe our approach to develop a simulation software application for the fully kinetic Vlasov equation which will be used to explore physics beyond the gyrokinetic model. Simulating the fully kinetic Vlasov equation requires efficient utilization of compute and storage capabilities due to the high dimensionality of the problem. In addition, the implementation needs to be extensibility regarding the physical model and flexible regarding the hardware for production runs. We start on the algorithmic background to simulate the 6-D Vlasov equation using a semi-Lagrangian algorithm. The performance portable software stack, which enables production runs on pure CPU as well as AMD or Nvidia GPU accelerated nodes, is presented. The extensibility of our implementation is guaranteed through the described software architecture of the main kernel, which achieves a memory bandwidth of almost 500 GB/s on a V100 Nvidia GPU and around 100 GB/s on an Intel Xeon Gold CPU using a single code base. We provide performance data on multiple node level architectures discussing utilized and further available hardware capabilities. Finally, the network communication bottleneck of 6-D grid based algorithms is quantified. A verification of physics beyond gyrokinetic theory for the example of ion Bernstein waves concludes the work.

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.

Continuum Vlasov simulations can be utilized for highly accurate modelling of fully kinetic plasmas. Great progress has been made recently regarding the applicability of the method in realistic plasma configurations. However, a reduction of the high computational cost that is inherent to fully kinetic simulations would be desirable, especially at high velocity space resolutions. For this purpose, low-rank approximations can be employed. The so far available low-rank solvers are restricted to either electrostatic systems or low dimensionality and can therefore not be applied to most space, astrophysical and fusion plasmas. In this paper we present a new parallel low-rank solver for the full six-dimensional electromagnetic Vlasov-Maxwell equations with a compression of the particle distribution function in velocity space. Special attention is paid to mass conservation and Gauss's law. The low-rank Vlasov solver is applied to standard benchmark problems of plasma turbulence and magnetic reconnection and compared to the full grid method. It yields accurate results at significantly reduced computational cost., Comment: 28 pages, 10 figures

Structure-preserving methods can be derived for the Vlasov-Maxwell system from a discretisation of the Poisson bracket with compatible finite-elements for the fields and a particle representation of the distribution function. These geometric electromagnetic particle-in-cell (GEMPIC) discretisations feature excellent conservation properties and long-time numerical stability. This paper extends the GEMPIC formulation in curvilinear coordinates to realistic boundary conditions. We build a de Rham sequence based on spline functions with clamped boundaries and apply perfect conductor boundary conditions for the fields and reflecting boundary conditions for the particles. The spatial semi-discretisation forms a discrete Poisson system. Time discretisation is either done by Hamiltonian splitting yielding a semi-explicit Gauss conserving scheme or by a discrete gradient scheme applied to a Poisson splitting yielding a semi-implicit energy-conserving scheme. Our system requires the inversion of the spline finite element mass matrices, which we precondition with the combination of a Jacobi preconditioner and the spectrum of the mass matrices on a periodic tensor product grid.

In this article we describe a unifying framework for variational electromagnetic particle schemes of spectral type, and we propose a novel spectral Particle-In-Cell (PIC) scheme that preserves a discrete Hamiltonian structure. Our work is based on a new abstract variational derivation of particle schemes which builds on a de Rham complex where Low's Lagrangian is discretized using a particle approximation of the distribution function. In this framework, which extends the recent Finite Element based Geometric Electromagnetic PIC (GEMPIC) method to a variety of field solvers, the discretization of the electromagnetic potentials and fields is represented by a de Rham sequence of compatible spaces, and the particle-field coupling procedure is described by approximation operators that commute with the differential operators in the sequence. In particular, for spectral Maxwell solvers the choice of truncated $L^2$ projections using continuous Fourier transform coefficients for the commuting approximation operators yields the gridless Particle-in-Fourier method, whereas spectral Particle-in-Cell methods are obtained by using discrete Fourier transform coefficients computed from a grid. By introducing a new sequence of spectral pseudo-differential approximation operators, we then obtain a novel variational spectral PIC method with discrete Hamiltonian structure that we call Fourier-GEMPIC. Fully discrete schemes are then derived using a Hamiltonian splitting procedure, leading to explicit time steps that preserve the Gauss laws and the discrete Poisson bracket associated with the Hamiltonian structure. These explicit steps share many similarities with standard spectral PIC methods. As arbitrary filters are allowed in our framework, we also discuss aliasing errors and study a natural back-filtering procedure to mitigate the damping caused by anti-aliasing smoothing particle shapes., Comment: 40 pages, 34 figures

In this article we apply a discrete action principle for the Vlasov--Maxwell equations in a structure-preserving particle-field discretization framework. In this framework the finite-dimensional electromagnetic potentials and fields are represented in a discrete de Rham sequence involving general finite element spaces, and the particle-field coupling is represented by a set of projection operators that commute with the differential operators. With a minimal number of assumptions which allow for a variety of finite elements and shape functions for the particles, we show that the resulting variational scheme has a general discrete Poisson structure and thus leads to a semi-discrete Hamiltonian system. By introducing discrete interior products we derive a second type of space discretization which is momentum preserving, based on the same finite elements and shape functions. We illustrate our method by applying it to spline finite elements, and to a new spectral discretization where the particle-field coupling relies on discrete Fourier transforms., Comment: 37 pages, 18 figures

This paper investigates subcycling of particle orbits in variational, geometric particle-in-cell methods addressing the Vlasov--Maxwell system in magnetized plasmas. The purpose of subcycling is to allow different time steps for different particle species and, ideally, time steps longer than the electron gyroperiod for the global field solves while sampling the local cyclotron orbits accurately. The considered algorithms retain the electromagnetic gauge invariance of the discrete action, guaranteeing a local charge conservation law, while the variational approach provides a bounded long-time energy behavior., Comment: 17 pages, 7 figures, accepted for publication in Physics of Plasmas

Numerical schemes that preserve the structure of the kinetic equations can provide stable simulation results over a long time. An electromagnetic particle-in-cell solver for the Vlasov-Maxwell equations that preserves at the discrete level the non-canonical Hamiltonian structure of the Vlasov-Maxwell equations has been presented in [Kraus et al. 2017]. Whereas the original formulation has been obtained for Cartesian coordinates, we extend the formulation to curvilinear coordinates in this paper. For the discretisation in time, we discuss several (semi-)implicit methods either based on a Hamiltonian splitting or a discrete gradient method combined with an antisymmetric splitting of the Poisson matrix and discuss their conservation properties and computational efficiency.

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 figures

This paper reports on an in-depth evaluation of the performance portability frameworks Kokkos and RAJA with respect to their suitability for the implementation of complex particle-in-cell (PIC) simulation codes, extending previous studies based on codes from other domains. At the example of a particle-in-cell model, we implemented the hotspot of the code in C++ and parallelized it using OpenMP, OpenACC, CUDA, Kokkos, and RAJA, targeting multi-core (CPU) and graphics (GPU) processors. Both, Kokkos and RAJA appear mature, are usable for complex codes, and keep their promise to provide performance portability across different architectures. Comparing the obtainable performance on state-of-the art hardware, but also considering aspects such as code complexity, feature availability, and overall productivity, we finally draw the conclusion that the Kokkos framework would be suited best to tackle the massively parallel implementation of the full PIC model.

This paper discusses energy-conserving time-discretizations for finite element particle-in-cell discretizations of the Vlasov--Maxwell system. A geometric spatially discrete system can be obtained using a standard particle-in-cell discretization of the particle distribution and compatible finite element spaces for the fields to discretize the Poisson bracket of the Vlasov--Maxwell model (see Kraus et al., J Plasma Phys 83, 2017). In this paper, we derive energy-conserving time-discretizations based on the discrete gradient method applied to an antisymmetric splitting of the Poisson matrix. Firstly, we propose a semi-implicit method based on the average-vector-field discretization of the subsystems. Moreover, we devise an alternative discrete gradient that yields a time discretization that can additionally conserve Gauss' law. Finally, we explain how substepping for fast species dynamics can be incorporated.

This work proposes a basis for improved throughput of matrix-free evaluation of discontinuous Galerkin symmetric interior penalty discretizations on hexahedral elements. The basis relies on ideas of Hermite polynomials. It is used in a fully discontinuous setting not for higher order continuity but to minimize the effective stencil width, namely to limit the neighbor access of an element to one data point for the function value and one for the derivative. The basis is extended to higher orders with nodal contributions derived from roots of Jacobi polynomials and extended to multiple dimensions with tensor products, which enable the use of sum factorization. The beneficial effect of the reduced data access on modern processors is shown. Furthermore, the viability of the basis in the context of multigrid solvers is analyzed. While a plain point-Jacobi approach is less efficient than with the best nodal polynomials, a basis change via sum-factorization techniques enables the combination of the fast matrix-vector products with effective multigrid constituents. The basis change is essentially for free on modern hardware because these computations can be hidden behind the cost of the data access., Comment: 20 pages, 8 figures

Gaussian kernels can be an efficient and accurate tool for multivariate interpolation. In practice, high accuracies are often achieved in the flat limit where the interpolation matrix becomes increasingly ill-conditioned. Stable evaluation algorithms for isotropic Gaussians (Gaussian radial basis functions) have been proposed based on a Chebyshev expansion of the Gaussians by Fornberg, Larsson & Flyer and based on a Mercer expansion with Hermite polynomials by Fasshauer & McCourt. In this paper, we propose a new stabilization algorithm for the multivariate interpolation with isotropic or anisotropic Gaussians derived from the generating function of the Hermite polynomials. We also derive and analyse a new analytic cut-off criterion for the generating function expansion that allows to automatically adjust the number of the stabilizing basis functions.

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.

This work presents algorithms for the efficient implementation of discontinuous Galerkin methods with explicit time stepping for acoustic wave propagation on unstructured meshes of quadrilaterals or hexahedra. A crucial step towards efficiency is to evaluate operators in a matrix-free way with sum-factorization kernels. The method allows for general curved geometries and variable coefficients. Temporal discretization is carried out by low-storage explicit Runge-Kutta schemes and the arbitrary derivative (ADER) method. For ADER, we propose a flexible basis change approach that combines cheap face integrals with cell evaluation using collocated nodes and quadrature points. Additionally, a degree reduction for the optimized cell evaluation is presented to decrease the computational cost when evaluating higher order spatial derivatives as required in ADER time stepping. We analyze and compare the performance of state-of-the-art Runge-Kutta schemes and ADER time stepping with the proposed optimizations. ADER involves fewer operations and additionally reaches higher throughput by higher arithmetic intensities and hence decreases the required computational time significantly. Comparison of Runge-Kutta and ADER at their respective CFL stability limit renders ADER especially beneficial for higher orders when the Butcher barrier implies an overproportional amount of stages. Moreover, vector updates in explicit Runge--Kutta schemes are shown to take a substantial amount of the computational time due to their memory intensity.

We present an algorithmic framework for matrix-free evaluation of discontinuous Galerkin finite element operators based on sum factorization on quadrilateral and hexahedral meshes. We identify a set of kernels for fast quadrature on cells and faces targeting a wide class of weak forms originating from linear and nonlinear partial differential equations. Different algorithms and data structures for the implementation of operator evaluation are compared in an in-depth performance analysis. The sum factorization kernels are optimized by vectorization over several cells and faces and an even-odd decomposition of the one-dimensional compute kernels. In isolation our implementation then reaches up to 60\% of arithmetic peak on Intel Haswell and Broadwell processors and up to 50\% of arithmetic peak on Intel Knights Landing. The full operator evaluation reaches only about half that throughput due to memory bandwidth limitations from loading the input and output vectors, MPI ghost exchange, as well as handling variable coefficients and the geometry. Our performance analysis shows that the results are often within 10\% of the available memory bandwidth for the proposed implementation, with the exception of the Cartesian mesh case where the cost of gather operations and MPI communication are more substantial.

Gaussian radial basis functions can be an accurate basis for multivariate interpolation. In practise, high accuracies are often achieved in the flat limit where the interpolation matrix becomes increasingly ill-conditioned. Stable evaluation algorithms have been proposed by Fornberg, Larsson & Flyer based on a Chebyshev expansion of the Gaussian basis and by Fasshauer & McCourt based on a Mercer expansion with Hermite polynomials. In this paper, we propose another stabilization algorithm based on Hermite polynomials but derived from the generating function of Hermite polynomials. The new expansion does not require a complicated choice of parameters and offers a simple extension to high-dimensional tensor grids as well as a generalization for anisotropic multivariate basis functions using the Hagedorn generating function.

Error estimates for the numerical solution of the master equation are presented. Estimates are based on adjoint methods. We find that a good estimate can often be computed without spending computational effort on a dual problem. Estimates are applicable to both settings with time-independent, and time-dependent propensity functions. The Finite State Projection algorithm reduces the dimensionality of the problem and time propagation is based on an Arnoldi exponential integrator, which in the time-dependent setting is combined with a Magnus method. Local error estimates are devised for the truncation of both the Magnus expansion and the Krylov subspace in the Arnoldi algorithm. An issue with existing methods is that error estimates for truncation of the state space depend on measuring a loss of probability mass in a way that is not usually compatible with the approximation of the exponential. We suggest an alternative error estimate that is compatible with a Krylov approximation to the matrix exponential. Finally, we apply the new error estimates to develop an adaptive simulation algorithm. Numerical examples demonstrate the benefits of the approach.

We present a novel framework for Finite Element Particle-in-Cell methods based on the discretization of the underlying Hamiltonian structure of the Vlasov-Maxwell system. We derive a semi-discrete Poisson bracket, which retains the defining properties of a bracket, anti-symmetry and the Jacobi identity, as well as conservation of its Casimir invariants, implying that the semi-discrete system is still a Hamiltonian system. In order to obtain a fully discrete Poisson integrator, the semi-discrete bracket is used in conjunction with Hamiltonian splitting methods for integration in time. Techniques from Finite Element Exterior Calculus ensure conservation of the divergence of the magnetic field and Gauss' law as well as stability of the field solver. The resulting methods are gauge invariant, feature exact charge conservation and show excellent long-time energy and momentum behaviour. Due to the generality of our framework, these conservation properties are guaranteed independently of a particular choice of the Finite Element basis, as long as the corresponding Finite Element spaces satisfy certain compatibility conditions., Comment: 57 Pages

In this article, we derive a semi-Lagrangian scheme for the solution of the Vlasov equation represented as a low-parametric tensor. Grid-based methods for the Vlasov equation have been shown to give accurate results but their use has mostly been limited to simulations in two dimensional phase space due to extensive memory requirements in higher dimensions. Compression of the solution via high-order singular value decomposition can help in reducing the storage requirements and the tensor train (TT) format provides efficient basic linear algebra routines for low-rank representations of tensors. In this paper, we develop interpolation formulas for a semi-Lagrangian solver in TT format. In order to efficiently implement the method, we propose a compression of the matrix representing the interpolation step and an efficient implementation of the Hadamard product. We show numerical simulations for standard test cases in two, four and six dimensional phase space. Depending on the test case, the memory requirements reduce by a factor $10^2-10^3$ in four and a factor $10^5-10^6$ in six dimensions compared to the full-grid method.

In this paper, we present a block-oriented scheme for adaptive mesh refinement based on summation-by-parts (SBP) finite difference methods and simultaneous-approximation-term (SAT) interface treatment. Since the order of accuracy at SBP-SAT grid interfaces is lower compared to that of the interior stencils, we strive at using the interior stencils across block-boundaries whenever possible. We devise a stable treatment of SBP-FD junction points, i.e. points where interfaces with different boundary treatment meet. This leads to stable discretizations for more flexible grid configurations within the SBP-SAT framework, with a reduced number of SBP-SAT interfaces. Both first and second derivatives are considered in the analysis. Even though the stencil order is locally reduced close to numerical interfaces and corner points, numerical simulations show that the locally reduced accuracy does not severely reduce the accuracy of the time propagated numerical solution. Moreover, we explain how to organize the grid and how to automatically adapt the mesh, aiming at problems of many variables. Examples of adaptive grids are demonstrated for the simulation of the time-dependent Schr\"odinger equation and for the advection equation.