Particle-in-cell methods (PIC)

Basics

PIC methods are well-suited for the discretization of conservation laws of the form

(1)\[\begin{split}\partial_t f + \nabla_q \cdot (\mathbf u(t, q) \, f) &= S(f) \qquad q \in \Omega \subset \mathbb R^n\,, \\[2mm] f(t=0) &= f_\textnormal{in}\,,\end{split}\]

where \(\Omega\subset \mathbb R^n\) is open, \(q\) denotes a set of suitable coordinates in \(\Omega\), and we assume appropriate boundary conditions on \(\partial\Omega\). Here, \(f: (t,q) \mapsto \mathbb R^+\) is a positive density or “volume form”, \(\mathbf u: (t,q) \mapsto \mathbb R^n\) is the vector-field defining the flow, and \(S(f)\) is a (nonlinear) source term. Equation (1) can be re-written as

\[\partial_t f + \mathbf u \cdot \nabla_q f = S(f) - f \nabla_q \cdot \mathbf u \,,\]

which represents an advection equation with a source term. The PIC ansatz for solving equation (1) is

(2)\[f \approx f_h(t, q) = \frac 1N \sum_{k=1}^N w_k(t)\, \delta(q - q_k(t))\,,\]

where \(N \gg 1\) is a large number and \(\delta(q)\) stands for the Dirac delta-distribution. Moreover, \(q_k: \mathbb R \to \Omega\) denotes the trajectory of marker \(k\) in the domain \(\Omega\) and \(w_k: \mathbb R \to \mathbb R\) stands for its (time-dependent) “weight”. The equations for \(q_k(t)\) and \(w_k(t)\) are derived from the following principle:

(3)\[\frac{\textrm{d}}{\textrm{d} t} \int_\Omega (f - f_h)\, \phi\,\textrm d q = O(N^{-1/2}) \qquad \forall \ \phi \in C^\infty(\Omega)\,.\]

The time derivative of the \(f_h\)-integral reads

(4)\[\begin{split}\frac{\textrm{d}}{\textrm{d} t} \int_\Omega f_h \phi\,\textrm d q &= \frac{\textrm{d}}{\textrm{d} t} \frac 1N \sum_{k=1}^N w_k\, \phi(q_k) \\[2mm] &= \frac 1N \sum_{k=1}^N \left[ \dot w_k\, \phi(q_k) + w_k\, \dot q_k \cdot \nabla_q \phi(q_k) \right]\,.\end{split}\]

Assuming that boundary terms vanish, the time derivative of the \(f\)-integral can be written as

\[\begin{split}\frac{\textrm{d}}{\textrm{d} t} \int_\Omega f \phi\,\textrm d q &= - \int_\Omega \nabla_q \cdot (\mathbf u \, f) \phi\,\textrm d q + \int_\Omega S(f) \phi\,\textrm d q \\[2mm] &= \int_\Omega f\,\mathbf u \cdot \nabla_q \phi\,\textrm d q + \int_\Omega S(f) \phi\,\textrm d q \,.\end{split}\]

These integrals are now viewed as Monte-Carlo integrals in the following way: assume \(s: (t,q) \mapsto \mathbb R^+\) to be a time-dependent probability distribution function (PDF) on \(\Omega\), and let \(\mathbb E(X)_s\) denote the expectation value of a random variable \(X\) distributed according to \(s\). Then,

\[\begin{split}\int_\Omega f\,\mathbf u \cdot \nabla_q \phi\,\textrm d q &= \int_\Omega \frac fs\,\mathbf u \cdot \nabla_q \phi\,s\,\textrm d q = \mathbb E\left(\frac fs\,\mathbf u \cdot \nabla_q \phi \right)_s\,, \\[2mm] \int_\Omega S(f) \phi\,\textrm d q &= \int_\Omega \frac 1s S(f) \phi\,s \,\textrm d q = \mathbb E\left(\frac 1s S(f) \phi \right)_s\,.\end{split}\]

At each time \(t\) we can estimate the above expectation values from the \(N\) samples \(q_k\):

(5)\[\begin{split}\mathbb E\left(\frac fs\,\mathbf u \cdot \nabla_q \phi \right)_s &\approx \frac 1N \sum_{k=1}^N \frac{f(q_k)}{s(q_k)}\,\mathbf u(t, q_k) \cdot \nabla_q \phi(q_k) + O(N^{-1/2})\,, \\[2mm] \mathbb E\left(\frac 1s S(f) \phi \right)_s &\approx \frac 1N \sum_{k=1}^N\frac{1}{s(q_k)}\,S(f(q_k)) \phi(q_k) + O(N^{-1/2})\,.\end{split}\]

We can compare this to equation (4) and since equation (3) holds for any \(\phi\), by comparing coefficients we deduce

\[\begin{split}w_k\, \dot q_k &= \frac{f(q_k)}{s(q_k)}\,\mathbf u(t, q_k)\,, \\[2mm] \dot w_k &= \frac{1}{s(q_k)}\,S(f(q_k))\,.\end{split}\]

Let us define the weights to be

(6)\[w_k := \frac{f(q_k)}{s(q_k)} \,,\]

which implies

(7)\[\begin{split}\dot q_k &= \mathbf u(t, q_k)\,, \\[2mm] \dot w_k &= \frac{1}{s(q_k)}\,S(f(q_k))\,.\end{split}\]

and also

(8)\[\dot w_k = \frac{\dot f(q_k)}{s(q_k)} - \frac{f(q_k)}{s(q_k)^2} \dot s(q_k) \,.\]

By comparing the equations (7) and (8) we find the equation that must be satisfied by the PDF \(s\):

\[\dot s(q_k) = \frac{s(q_k)}{f(q_k)} \Big[ \dot f(q_k) - S(f(q_k)) \Big]\,.\]

Since, from equation (1),

\[\dot f(q_k) = S(f(q_k)) - f(q_k) \nabla_q \cdot \mathbf u(t, q_k)\,,\]

we deduce

\[\dot s(q_k) = - s(q_k) \nabla_q \cdot \mathbf u(t, q_k) \,.\]

Therefore, with the choice of the weights made in equation (6), the PDF \(s\) must satisfy

\[\partial_t s + \nabla_q \cdot (\mathbf u(t, q) \, s) = 0\,.\]

This means that \(s\) is a density (or volume form) advected by the flow of \(\mathbf u\). For incompressible flow, \(\nabla_q \cdot \mathbf u = 0\), the PDF \(s\) is constant along the flow of \(\mathbf u\) and the PIC solution (2) is obtained by solving

(9)\[\begin{split}\dot q_k &= \mathbf u(t, q_k) \qquad &&q_k(0) = q_{k0}\,, \\[2mm] \dot w_k &= \frac{1}{s(t=0, q_{k0})}\,S(f) \qquad &&w_k(0) = \frac{f(t=0, q_{k0})}{s(t=0, q_{k0})}\,,\end{split}\]

where the \(q_{k0}\) are drawn according to the PDF \(s(t=0)\). On the other hand, for compressible flow \(\nabla_q \cdot \mathbf u \neq 0\) the weights update is more complicated:

(10)\[\begin{split}\dot q_k &= \mathbf u(t, q_k) \qquad &&q_k(0) = q_{k0}\,, \\[2mm] \dot w_k &= \frac{w_k}{f(t, q_k)}\,S(f) \qquad &&w_k(0) = \frac{f(t=0, q_{k0})}{s(t=0, q_{k0})}\,.\end{split}\]

Note that when there is no source, \(S(f) = 0\), the weights are constant for compressible and incompressible flow.

PIC methods are popular for solving these types of equations because

1. they do not suffer from the curse of dimensionality for large \(n\)

2. they and can be easily parallelized.

However, a disadvantage is their slow convergence \(\sim 1/\sqrt N\) with the number of particles (also called markers) and the associated noise that is present in simulations. Often times noise-reduction techniques such as the Control variate method have to be employed to reduce simulation cost.

Monte-Carlo integrals

At the core of the PIC method is the approximation of integrals

\[I = \int_{\Omega} f(q)\, A(q) \,\textnormal d^n q\,.\]

by Mont-Carlo (MC) integration. The first expecation value in equation (5) is an example of this. In Struphy, such integrals are calculated in non-Cartesian coordinates \(q \in \mathbb R^n\), and we denote by

\[\textnormal d^n q = |J| \textnormal d q\,,\]

the volume element, including the Jacobian determinant from the mapping to Cartesian (straight) space. Thus, we could also write the integral as

(11)\[I = \int_{\Omega} f(q)\, A(q) |J|\,\textnormal d q\,,\]

We have seen that \(s\) satisfies the transport equation s_eq, and that it is normalized, thus

\[\int_{\Omega} s |J|\,\textnormal d q = 1 \qquad \forall t\,, \qquad \quad s >0\,.\]

In MC integration, we view the integral \(I\) as an expectation value:

(12)\[I = \int_{\Omega} \frac{f(q)}{s(q)}\, A(q)\, s(q) |J|\,\textnormal d q = \mathbb E\left( \frac{f}{s} A\right)_s\,,\]

with respect to the PDF \(s|J|\). Having access to \(N\) samples \(q_k(t)\) of \(s\) at each time \(t\) allows us to estimate this expectation via

(13)\[I = \mathbb E\left( \frac{f}{s} A\right)_s \approx \frac 1N \sum_{k=1}^N \frac{f(q_k)}{s(q_k)} A(q_k) = \frac 1N \sum_{k=1}^N w_k A(q_k)\,,\]

using the definition of the weights in equation (6). The error in this approximation is of order \(N^{-1/2}\).

The weights are implemented in struphy.pic.base.Particles.initialize_weights(). Equation (13) can be obtained directly from (12) by inserting the PIC ansatz (2), with \(f|J| \approx f_h\).

In Struphy, Monte-Carlo integrals of the form (13) are implemented via accums.

PIC algorithm

The PIC algorithm can be summarized in the following steps:

1. Draw \(N\gg 1\) markers \((q_{k0})_k\) according to \(s_{\textnormal{in}}|J|\). The initial marker distribution is implemented in struphy.pic.base.Particles.draw_markers(). In Struphy, it is always a Gaussian in velocity space (see docstring).

2. Compute the weights \((w_{k0})_k\) according to equation (6). The initial condition of the kinetic equation enters here.

3. Solve the characteristic equations (9) for each marker.

4. Whenever necessary, compute integrals according to (13).

Particle binning

The aim of particle binning is to obtain a “grid-representation” of either \(f^0(t, q)\) or \(f^\textrm{vol}(t, q)\), which are represented by markers in Struphy. The approximation is obtained by integration of \(f^{0/\textrm{vol}}(t, q)\) over small volumes of phase space, called “bins”. These bins are not necessarily related to the FE grid used for fluid/field variables (but they can be). Moreover, bins do not have to be defined in all \(n\) directions of the phase space; one could do just 1D or 2D binning, integrating over the other coordinates entirely.

For example, let’s say we want to do 2D binning in the \((\eta_1,v_1)\) subspace of the full \((\eta,v) \in \mathbb R^6\) phase space. In this case we want an approximate representation of the “reduced” 2D density

\[f^{k,\textnormal{red}}(t, \eta_1, v_1) := \int_{(0, 1)^2}\int_{\mathbb R^2} f^k(t, \eta, v) \,\textnormal d \eta_2 \textnormal d \eta_3 \textnormal d v_2 \textnormal dv_3\,,\]

where \(k \in \{0, \textrm{vol}\}\) such that \(f^k\) is either the 0-form or the volume form in logical space (here integrated over a submanifold, which is technically not allowed - but serves to get an approximation via numerical binning). If we define bin edges \((\eta_{1i})_i\) and \((v_{1j})_j\) in the \((\eta_1, v_1)\) subspace, we can approximate \(f^{k,\textnormal{red}}\) restricted to a bin \(\Omega_{ij} = (\eta_{1i}, \eta_{1i+1}) \times (v_{1j}, v_{1j+1})\) by its mean value over over that bin:

(14)\[f^{k,\textnormal{red}}(t, \eta, v)\Big|_{\Omega_{ij}} \approx R^k_{ij}(t) := \frac{1}{|\Omega_{ij}|} \int_{\Omega_{ij}} f^{k,\textnormal{red}}(t, \eta_1, v_1)\,\textnormal d\eta_1 \textnormal d v_1\,.\]

where \(|\Omega_{ij}|\) is the volume of the \(ij\)-th bin (in 2D):

\[|\Omega_{ij}| = \int_{\Omega_{ij}} \textnormal d\eta_1 \textnormal d v_1 = (\eta_{1i+1} - \eta_{1i})(v_{1j+1} - v_{1j})\,.\]

Note that \(R^k_{ij}\) has the same dimension as \(f^{k,\textnormal{red}}\), due to the one-over-volume prefactor. The integrals \(R^k_{ij}\) are then approximated by Monte-Carlo integration, using the indicator function \(\textnormal{id}_{ij}(x_1, v_1)\) of the \(ij\)-th bin (it is 1 inside \(\Omega_{ij}\) and zero otherwise). Depending on the value of \(k \in \{0, \textrm{vol}\}\), there are two different integrals to approximate:

\begin{align} R^0_{ij}(t) &= \frac{1}{|\Omega_{ij}|}\int_{(0, 1)^3} \int_{\mathbb R^3} \textnormal{id}_{ij}(x_1, v_1)\, f^0(t, \eta, v)\,\textnormal d\eta\, \textnormal d v \\ &= \frac{1}{|\Omega_{ij}|} \int_{(0, 1)^3} \int_{\mathbb R^3} \textnormal{id}_{ij}(x_1, v_1)\, \frac{f^0}{s^0 |\textnormal{det} J_F|}\,s^\textrm{vol}\,\textnormal d\eta\, \textnormal d v \\ &\approx \frac{1}{N|\Omega_{ij}|} \sum_{k \in \Omega_{ij}} w_{k} \frac{1}{|\textnormal{det} J_F|(\eta_k(t), v_k(t))}\,, \\[2mm] R^\textrm{vol}_{ij}(t) &= \frac{1}{|\Omega_{ij}|}\int_{(0, 1)^3} \int_{\mathbb R^3} \textnormal{id}_{ij}(x_1, v_1)\, f^\textrm{vol}(t, \eta, v)\,\textnormal d\eta\, \textnormal d v \\ &= \frac{1}{|\Omega_{ij}|} \int_{(0, 1)^3} \int_{\mathbb R^3} \textnormal{id}_{ij}(x_1, v_1)\, \frac{f^0}{s^0}\,s^\textrm{vol}\,\textnormal d\eta\, \textnormal d v \\ &\approx \frac{1}{N|\Omega_{ij}|} \sum_{k \in \Omega_{ij}} w_{k} \,. \end{align}

Here, the sum runs over all markers that are in the bin \(\Omega_{ij}\). The values \(R^0_{ij}\) computed in (15) are discrete approximations of \(f^0\), and values \(R^\textrm{vol}_{ij}\) are discrete approximations of \(f^\textrm{vol}\). This binning algorithm is implemented in struphy.pic.base.Particles.binning().

Control variate method

This noise reduction technique relies on

1. The ability to compute (without Monte-Carlo) the integral (11) for a known distribution function \(f^0 = \mathcal M^0\):

\[I_{\mathcal M} := \int_{\mathbb R^n} \mathcal M^0 (t,q)\, A^0(t,q) |J_F(q)|\,\textnormal d q = \int_{\mathbb R^n} \mathcal M^n (t,q)\, A^0(t,q) \,\textnormal d q\,.\]

2. The assumption that \(f^0(t)\) is close to \(\mathcal M^0(t)\) for all \(t\).

In that case we can decompose (11) as

\begin{align} I &= \int_{\mathbb R^n} f^0 \, A^0 |J_F|\,\textnormal d q\,, \\[1mm] &= \int_{\mathbb R^n} (f^0 - \mathcal M^0) \, A^0 |J_F|\,\textnormal d q + I_{\mathcal M}\,, \\[1mm] &= \int_{\mathbb R^n} \frac{(f^0 - \mathcal M^0)}{s^0} \, A^0 s^\textrm{vol}\,\textnormal d q + I_{\mathcal M}\,, \\[1mm] &\approx \frac{1}{N} \sum_{k=1}^N w_{k}(t)\, A^0 (t, q_k(t)) + I_{\mathcal M}\,, \end{align}

where we have introduced the time-dependent weights

\[w_k(t) = w_{k0} - \frac{\mathcal M^0(t, q_k(t))}{s^0_{\textnormal{in}}(q_{k0})}\,.\]

This algorithm is implemented in struphy.pic.base.Particles.update_weights().