GitHub Repository for this project
PIC-Simulation
Nigel-Shen ā¢ Updated Sep 2, 2022
šŖ Backgroundš±Ā Particle-in-cellšĀ Particle-in-FourieršĀ Free space Poisson solverš”Ā Our Work: combining PIF with the free space solver
šŖ Background
Traditionally, when we describe the fluid motion we use the Eulerās equation
which is essentially a transport equation in 2D / 3D, under the assumption that the fluid is fully collisional, meaning that locally the change of particle velocities is smooth.
āMacroscopically small, microscopically largeāā¦
This is not the case when we encounter plasma, a very heated and ionized state of matter. In this case, a particle might have neighbors whose velocities are drastically different ā in other words, in a spatial location we donāt have one specific velocity value, but a distribution.
In order to describe the dynamics of plasma, we need a PDE that describes the distribution of plasma not only in physical space but also in velocity space. This is the so called Vlasovās equation (Boltzmannās equation if youāre a physicist):
where is the distribution of plasma in the phase space. The derivation of Vlasovās equation is similar to that of Eulerās ā just start with , and use the chain rule.
The charge density is defined as:
where is the total charge of the system. A few conservation laws include:
where
One of the problems in numerically solving this equation is that it lives in a space up to 6 dimensions, which significantly slows down the numerical algorithms (curse of dimensionality) if we simply use tools like finite difference or spectral methods.
š±Ā Particle-in-cell
Instead, Particle-In-Cell (PIC) Method solves the equation by introducing finite-size particles, each has its own charge, mass, location (smoothed by a shape function), and velocity. We can then simulate the motion of these particles by simply calculating the electric field ( complexity), and do a particle push. This reduces our model to a 3-dimensional problem.
In order to solve for the electric field in PIC, we spread particle charges onto gridpoints using the shape function (similar to Immersed Boundary), solve the Poissonās equation on the grid, and interpolate the values back to the particle positions:
However, this is not a very elegant way, because in this interpolation scheme, the charges on the grid points are a discretized version of the true particle distribution, and vice versa. Some information was lost during this process. This causes the energy to always increase over the time of simulation, making the entire scheme unreliable for long time simulations.
šĀ Particle-in-Fourier
Recently, a new scheme called Particle-in-Fourier has been proposed. It utilizes nonuniform FFT to interpolate particles directly into the Fourier space:
The field we obtain in Fourier space is spectral accurate, and it acts on the particles directly, with no information lost or truncated during interpolation. As a numerical algorithm, however, there must be some truncation because computers only work with discrete numbers. It turns out the shape function is Fourier truncated to become a band-limited function because of finite number of Fourier modes (shape function can be exact in PIC), but this does not matter so much in our implementation.
One can prove that unlike PIC, PIF provides exact energy conservation as . This has been checked by doing some analytical calculation. We then make an educated guess that the convergence of energy should be of the same order as our time integrator. (In our work, we proved second order convergence under leapfrog)
Another interpretation: For a band-limited function , the solution we get from the spectral method is not spectral accurate, but exactly correct (assuming we have enough number of modes). Thatās why the field solver does not introduce any numerical error apart from machine error, so all conservation properties should preserve exactly.
šĀ Free space Poisson solver
Free space Poisson solvers are difficult to construct, as the solution is always entire in and obviously computers cannot deal with infinity. Often people use precomputed Greenās functions which are compact supported in the space and convolute it with the source .
Suppose the Greenās function in free space is , then in our computational domain , it makes no difference convoluting with or with where
is a radial symmetric, compact supported function, and . It turns out the Fourier transform of can be computed analytically and decays exponentially:
We can thus use to multiply with to get the free space solution of the Poissonās equation. Moreover, this method is spectral accurate.
It seems a coincidence that this method is spectral accurate. We know that the Greenās function is very singular, and it seems that the mollified function brings even more discontinuity, which seems hard to be handled using FFT. However, analytical calculation shows that the Fourier coefficients of decays exponentially, which is amazing.
š”Ā Our Work: combining PIF with the free space solver
Observe that both PIF and the free space Poisson solver work directly in the Fourier space. It is thus natural to combine these two schemes to create a spectral accurate free space PIF solver.
There is some math involved, but the heuristics is to properly define these two convolutional kernels:
This involves a more careful choice of , and changing the order of some computation steps by taking advantage of the linearity of the convolution and differentiation operators.
Moreover, we can generalize this to arbitrary Dirichlet boundary conditions (and maybe Neumann, but we didnāt implement this) by applying an additional Laplace solver, which can be done using immersed boundary / interface, or boundary integral method. In our case, we donāt care too much about the convergence near boundary, since most particles stays at the center of the computational domain, so we use boundary integral, which provides spectral accuracy.
Thereās some subtlety which involves convoluting of the homogeneous solution with shape function in a boundary integral method, which cannot be done in the usual manner. We took advantage of the symmetry of and the mean value property of the solution to maintain spectral accuracy (which I personally think is rather beautiful šš»). See our incoming preprint for more detailsā¦
