Abstract
Kinetic equations model distributions of particles in position-velocity phase space. Often, one is interested in studying the long-time behavior of particles in high-collisional regimes in which an approximate (advection)-diffusion model holds. In this paper we consider the diffusive scaling. Classical particle-based techniques suffer from a strict time-step restriction in this limit, to maintain stability. Asymptotic-preserving schemes avoid this problem, but introduce an additional time discretization error, possibly resulting in an unacceptably large bias for larger time steps. Here, we present and analyze a multilevel Monte Carlo scheme that reduces this bias by combining estimates using a hierarchy of different time step sizes. We demonstrate how to correlate trajectories from this scheme, using different time steps. We also present a strategy for selecting the levels in the multilevel scheme. Our approach significantly reduces the computation required to perform accurate simulations of the considered kinetic equations, compared to classical Monte Carlo approaches.
Similar content being viewed by others
References
Anderson, D.F., Higham, D.J.: Multilevel Monte Carlo for continuous time Markov chains, with applications in biochemical kinetics. Multiscale Model. Simul. 10(1), 146–179 (2012)
Bennoune, M., Lemou, M., Mieussens, L.: Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier–Stokes asymptotics. J. Comput. Phys. 227(8), 3781–3803 (2008)
Bhatnagar, P.L., Gross, E.P., Krook, M.: A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. 94(3), 511–525 (1954)
Birdsall, C.K., Langdon, A.B.: Plasma Physics via Computer Simulation. Series in Plasma Physics and Fluid Dynamics. Taylor & Francis, New York (2004)
Blondeel, P., Robbe, P., Van Hoorickx, C., François, S., Lombaert, G., Vandewalle, S.: p-refined multilevel quasi-Monte Carlo for Galerkin finite element methods with applications in civil engineering. Algorithms 13(5), 110 (2020)
Boscarino, S., Pareschi, L., Russo, G.: Implicit–explicit Runge–Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit. SIAM J. Sci. Comput. 35(1), A22–A51 (2013)
Buet, C., Cordier, S.: An asymptotic preserving scheme for hydrodynamics radiative transfer models. Numer. Math. 108(2), 199–221 (2007)
Cercignani, C.: The Boltzmann Equation and Its Applications, Applied Mathematical Sciences, vol. 67. Springer, New York (1988)
Cliffe, K.A., Giles, M.B., Scheichl, R., Teckentrup, A.L.: Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients. Comput. Vis. Sci. 14(1), 3–15 (2011)
Crestetto, A., Crouseilles, N., Lemou, M.: A particle micro–macro decomposition based numerical scheme for collisional kinetic equations in the diffusion scaling. Commun. Math. Sci. 16(4), 887–911 (2018)
Crouseilles, N., Lemou, M.: An asymptotic preserving scheme based on a micro–macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits. Kinet. Relat. Mod. 4(2), 441–477 (2011)
Degond, P., Dimarco, G., Pareschi, L.: The moment-guided Monte Carlo method. Int. J. Numer. Methods Fluids 67(2), 189–213 (2011)
Dimarco, G., Pareschi, L.: Hybrid multiscale methods II. Kinetic equations. Multiscale Model. Simul. 6(4), 1169–1197 (2008)
Dimarco, G., Pareschi, L.: Fluid solver independent hybrid methods for multiscale kinetic equations. SIAM J. Sci. Comput. 32(2), 603–634 (2010)
Dimarco, G., Pareschi, L.: High order asymptotic-preserving schemes for the Boltzmann equation. C. R. Math. 350(9–10), 481–486 (2012)
Dimarco, G., Pareschi, L.: Numerical methods for kinetic equations. Acta Numer. 23, 369–520 (2014)
Dimarco, G., Pareschi, L., Samaey, G.: Asymptotic-preserving Monte Carlo methods for transport equations in the diffusive limit. SIAM J. Sci. Comput. 40(1), A504–A528 (2018)
Giles, M.B.: Multilevel Monte Carlo path simulation. Oper. Res. 56(3), 607–617 (2008)
Giles, M.B.: Multilevel Monte Carlo methods. Acta Numer. 24, 259–328 (2015)
Goldstein, S.: On diffusion by discontinuous movements, and on the telegraph equation. Q. J. Mech. Appl. Math. 4(2), 129–156 (1951)
Gosse, L., Toscani, G.: An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations. C. R. Math. 334(4), 337–342 (2002)
Haji-Ali, A.L., Nobile, F., Tempone, R.: Multi-index Monte Carlo: when sparsity meets sampling. Numer. Math. 132(4), 767–806 (2016)
Heinrich, S.: Multilevel Monte Carlo methods. In: Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), vol. 2179, pp. 58–67. Springer (2001)
Hoel, H., Law, K.J.H., Tempone, R.: Multilevel ensemble Kalman filtering. SIAM J. Numer. Anal. 54(3), 1813–1839 (2016)
Jin, S.: Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations. SIAM J. Sci. Comput. 21(2), 441–454 (1999)
Jin, S., Pareschi, L., Toscani, G.: Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations. SIAM J. Numer. Anal. 35(6), 2405–2439 (1998)
Jin, S., Pareschi, L., Toscani, G.: Uniformly accurate diffusive relaxation schemes for multiscale transport equations. SIAM J. Numer. Anal. 38(3), 913–936 (2000)
Klar, A.: An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit. SIAM J. Numer. Anal. 35(3), 1073–1094 (1998)
Klar, A.: A numerical method for kinetic semiconductor equations in the drift-diffusion limit. SIAM J. Sci. Comput. 20(5), 1696–1712 (1999)
Lapeyre, B., Pardoux, É., Sentis, R., Craig, A.W., Craig, F.: Introduction to Monte Carlo Methods for Transport and Diffusion Equations, vol. 6. Oxford University Press, Oxford (2003)
Larsen, E.W., Keller, J.B.: Asymptotic solution of neutron transport problems for small mean free paths. J. Math. Phys. 15(1), 75–81 (1974)
Lemou, M., Mieussens, L.: A new asymptotic preserving scheme based on micro–macro formulation for linear kinetic equations in the diffusion limit. SIAM J. Sci. Comput. 31(1), 334–368 (2008)
Løvbak, E., Samaey, G., Vandewalle, S.: A multilevel Monte Carlo asymptotic-preserving particle method for kinetic equations in the diffusion limit. In: L’Ecuyer, P., Tuffin, B. (eds.) Monte Carlo and Quasi-Monte Carlo Methods 2018, pp. 383–402. Springer (2020)
Mortier, B., Baelmans, M., Samaey, G.: Kinetic-diffusion asymptotic-preserving Monte Carlo algorithms for plasma edge neutral simulation. Contrib. Plasma Phys. 60, e201900134 (2019)
Naldi, G., Pareschi, L.: Numerical schemes for hyperbolic systems of conservation laws with stiff diffusive relaxation. SIAM J. Numer. Anal. 37(4), 1246–1270 (2000)
Pareschi, L.: Hybrid multiscale methods for hyperbolic and kinetic problems. ESAIM Proc. 15, 87–120 (2005)
Pareschi, L., Caflisch, R.E.: An implicit Monte Carlo method for rarefied gas dynamics. J. Comput. Phys. 154(1), 90–116 (1999)
Pareschi, L., Russo, G.: An introduction to Monte Carlo method for the Boltzmann equation. ESAIM Proc. 10, 35–75 (2001)
Pareschi, L., Trazzi, S.: Numerical solution of the Boltzmann equation by time relaxed Monte Carlo (TRMC) methods. Int. J. Numer. Methods Fluids 48(9), 947–983 (2005)
Pope, S.B.: A Monte Carlo method for the PDF equations of turbulent reactive flow. Combust. Sci. Technol. 25(5–6), 159–174 (1981)
Rousset, M., Samaey, G.: Simulating individual-based models of bacterial chemotaxis with asymptotic variance reduction. Math. Models Methods Appl. Sci. 23(12), 2155–2191 (2011)
Van Barel, A., Vandewalle, S.: Robust optimization of PDEs with random coefficients using a multilevel Monte Carlo method. SIAM/ASA J. Uncertain. Quantif. 7(1), 174–202 (2019)
Acknowledgements
We thank the anonymous reviewer for their substantial suggestions which have significantly improved this article. Emil Løvbak is funded by the Research Foundation - Flanders (FWO) under SB-Fellowship number 1SB1919N. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government – department EWI.
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Conflict of interest
The authors declare that they have no conflict of interest.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
A preliminary version of this paper, cited as [33], was published in Monte Carlo and Quasi-Monte Carlo Methods 2018
Supplementary Information
Below is the link to the electronic supplementary material.
Appendices
Computing the covariance sums needed to prove Lemma 5
1.1 \(\displaystyle \sum _{m=0}^{M-1} {{\,\mathrm{\mathrm {Cov}}\,}}\left( \Delta T^{n,m}_{p,\ell }, \Delta T^{n}_{p,\ell -1} \right) \)
The covariance between the coarse increment n and a fine sub-increment (n, m), given Lemma 4, can be written out as
To calculate the expected value in (48), we need to consider the probabilities of coupled collisions taking place in the correlated simulations. If a collision takes place in a given set of M fine increments, the probability that the correlated coarse simulation time step will also simulate a collision is given by
From Lemma 1, we know that it is not possible for a collision to take place in the coarse simulation, without a collision taking place in the fine simulation. This leaves three possibilities, when considering collision behavior in coarse time step \(n-1\):
-
Both at level \(\ell -1\) and at level \(\ell \), no collision occurred in time step \(n-1\). In this case, time step \(n-1\) will not affect the correlation of the velocities between the simulations. If the velocities were correlated at the beginning of time step \(n-1\), they will still be so at the end of the time step, and vice versa. In this case, we thus need to look at step \(n-2\), and so on, until we reach a past time step that satisfies one of the following two cases.
-
A collision occurred at level \(\ell \) in time step \(n-1\), but not at level \(\ell -1\). In this case, a new \({\bar{V}}^{n-1,m,*}_{p,\ell }\) was drawn from \({{\mathcal {M}}}(v)\), for some m, independently of the value of \({\bar{V}}^{n-1,*}_{p,\ell -1}\). Because all sampled velocities are independent, there is no correlation between \({\bar{V}}^{n,m-1}_{p,\ell }\) and \({\bar{V}}^{n-1}_{p,\ell -1}\), making the expected value of their product zero by (10).
-
A collision occurred both at level \(\ell \) and at level \(\ell -1\) in time step \(n-1\). We know that \({\bar{V}}^{n,0}_{p,\ell } = {\bar{V}}^{n-1,M,*}_{p,\ell } = {\bar{V}}^{n}_{p,\ell -1} = {\bar{V}}^{n-1,*}_{p,\ell -1}\) by (29). For every other fine simulation sub-step \(m>1\) in time step n, there is a probability \(1-p_{nc,{\varDelta } t_\ell }^{m}\) that a collision has taken place in at least one of the steps (n, 0) through (n, m). These collisions will not affect the coarse simulation until time step \(n+1\), where \({\bar{V}}^{n,*}_{p,\ell }\) is used to generate a new velocity. So we have that \({\bar{V}}^{n,m}_{p,\ell } = {\bar{V}}^{n}_{p,\ell -1}\) with a probability \(p_{nc,{\varDelta } t_\ell }^{m}\), otherwise the two velocities are uncorrelated.
From the above list of possibilities, we conclude that the random variables \({\bar{V}}^{n,m}_{p,\ell }\) and \({\bar{V}}^{n}_{p,\ell -1}\) are equal if both of the following are true:
-
1.
The last simulation step at level \(\ell \) that underwent a collision is a sub-step of a simulation step \(n^\prime \) at level \(\ell -1\) which also underwent a collision.
-
2.
The coarse simulation step \(n^\prime \) is not the current coarse step n.
Otherwise, \({\bar{V}}^{n,m}_{p,\ell }\) and \({\bar{V}}^{n}_{p,\ell -1}\) are uncorrelated.
Making use of (10) and (49), we have that
where the final equality is a fortuitous coincidence that will shorten some expressions in Sect. 5.
1.2 \(\displaystyle \sum _{m=0}^{M-2} \sum _{m^\prime = m+1}^{M-1}{{\,\mathrm{\mathrm {Cov}}\,}}\left( \Delta T^{n,m}_{p,\ell }, \Delta T^{n,m^\prime }_{p,\ell } \right) \)
As in Sect. A.1, we start by writing out the covariance between subsequent fine sub-increments (n, m) and \((n,m^\prime )\) as
Calculating the sum of the covariance between subsequent fine increments is straightforward. To simplify notation, we introduce \(\Delta m\) as shorthand for \(m^\prime - m\). As the collision probability is constant across time steps, the covariance is given by
Making use of (10), we get that the r.h.s. of (53) is equal to
By making use of the identity
and noting that \(\frac{p_{nc,{\varDelta } t_\ell }}{p_{c,{\varDelta } t_\ell }}=\frac{\varepsilon ^2}{\Delta t_\ell }\), we work out (54) as
1.3 \(\displaystyle \sum _{n=0}^{N-2} \sum _{n^\prime =n+1}^{N-1}{{\,\mathrm{\mathrm {Cov}}\,}}\left( \Delta _{T,n} , \Delta _{T,n^\prime } \right) \)
Making use of Lemma 4, we write out the covariance of differences of transport increments for a given n and \(n^\prime \) as
We now calculate each term on the r.h.s. of (57) separately. We start by writing the expressions in terms of \(\Delta n\), as a shorthand for \(n^\prime - n\). In the remainder of the section we assume \(n^\prime > n\), without loss of generality.
The increment covariance at level \(\ell \) \((\text {I})\) is calculated from a similar starting point to (52):
We now substitute the double summation in (58) with a two summations over \(\mathrel {{\varDelta } m}= |m^\prime - m|\):
We now again make use of the identity (55), giving that (59) equals
Calculating the covariance of increments at level \(\ell -1\) (II) is also straightforward:
The expression for the covariance between increments at level \(\ell \) and level \(\ell -1\) at differing time steps n and \(n'\) depends on the relative position of the increments, i.e., whether the increment at level \(\ell \) comes before that at level \(\ell -1\), or not. If the fine increment at level \(\ell \) comes first (III), we need to calculate
To calculate the expected value of the r.h.s. of (62), we list the possible simulation behaviors at level \(\ell -1\) in time step n, relative to the fine sub-step m:
-
No collision occurs in the simulation at level \(\ell \). This situation occurs with probability \(p_{nc,{\varDelta } t_\ell }^M.\) In this case, we have already established that no collision occurs at level \(\ell -1\). The probability of the simulations being correlated at the end of time step n is therefore equal to the probability of them being correlated at the start of time step n, given by (49).
-
No collision occurs in sub-steps 0 through \(m-1\), but at least one collision occurs in sub-steps m through \(M-1\). This situation occurs with probability
$$\begin{aligned} p_{nc,{\varDelta } t_\ell }^{m} \left( 1 - p_{nc,{\varDelta } t_\ell }^{M-m} \right) = p_{nc,{\varDelta } t_\ell }^{m} - p_{nc,{\varDelta } t_\ell }^{M}. \end{aligned}$$If the simulation at level \(\ell -1\) also has a collision, then \({\bar{V}}^{n+1}_{p,\ell -1}\) and \({\bar{V}}^{n,m}_{p,\ell }\) are independent by (29). If no collision occurred in the simulation at level \(\ell -1\), which is the case with probability
$$\begin{aligned} \frac{p_{nc,{\varDelta } t_{\ell -1}}-p_{nc,{\varDelta } t_\ell }^M}{1-p_{nc,{\varDelta } t_\ell }^M}, \end{aligned}$$then \({\bar{V}}^n_{p,\ell -1} = {\bar{V}}^{n,m}_{p,\ell }\), if the trajectories were correlated at the beginning of time step n. This is the case with probability (49).
-
At least one collision occurs in sub-steps 0 through \(m-1\), but none occur in sub-steps m through M. This occurs with probability
$$\begin{aligned} \left( 1 - p_{nc,{\varDelta } t_\ell }^{m} \right) p_{nc,{\varDelta } t_\ell }^{M-m} = p_{nc,{\varDelta } t_\ell }^{M-m} - p_{nc,{\varDelta } t_\ell }^M. \end{aligned}$$By (29) we know that \({\bar{V}}^{n,*}_{p,\ell -1} = {\bar{V}}^{n,m}_{p,\ell }\), if a collision also occurs in time step n of the simulation at level \(\ell -1\). This happens with probability (49).
-
Collisions happen both before sub-step m and during or afterwards. In this case, there is no correlation between \({\bar{V}}^{n,*}_{p,\ell -1}\) and \({\bar{V}}^{n,m}_{p,\ell }\).
By adding the non-zero contributions from these four cases, we calculate the probability that \({\bar{V}}^{n,m}_{p,\ell }\) is correlated with \({\bar{V}}^{n,*}_{p,\ell -1}\). To get the probability of the correlation holding until \({\bar{V}}^{n^\prime }_{p,\ell -1}\), we multiply the sum by \(p_{nc,{\varDelta } t_{\ell -1}}^{\mathrel {{\varDelta } n}-1}\). Given the properties (10), we can state that the right hand side sum in (62) is given by
Making us of the fact that \(\frac{p_{c,{\varDelta } t_{\ell -1}}}{p_{nc,{\varDelta } t_{\ell -1}}}=\frac{\Delta t_{\ell -1}}{\varepsilon ^2}\) we can then write the r.h.s. of (63) as
To calculate the covariance of coarse time steps preceding fine time steps (IV), fewer calculations are needed. The two time steps are correlated if the trajectories were correlated at the start of coarse simulation time step n and this correlation was not lost due to a collision in the fine simulation between time steps (n, 0) through \((n^\prime ,m^\prime -1)\)
Making use of expressions (60) through (65), we get
\( \text {with} \qquad \sigma _1= \varepsilon ^2 {\tilde{v}}^2 \left( p_{nc,{\varDelta } t_\ell }^{1-M} + p_{nc,{\varDelta } t_\ell }^{M+1} - 2 p_{nc,{\varDelta } t_\ell }\right) - \Delta t_{\ell -1}^2 {\tilde{v}}_{\Delta t_{\ell -1}}^2 \qquad \text {and} \)
Theoretical background to Section 6.2
Based on Section 2.6 in [19] it can be shown that it is beneficial to leave out level 0 if
with \(C_\ell \) the cost of a sample at level \(\ell \). Our analytical results from Sect. 6.1 can only strictly be applied in (67) if the quantity of interest is the expected particle position or velocity, as applying Lipschitz constants gives independent upper bounds for the positive and negative terms, hence we will also back up our claims with numerical experiments. Here, we only consider the particle position variance, as \({{\,\mathrm{{\mathbb {V}}}\,}}\left[ X_\ell -X_{\ell -1}\right] \gg {{\,\mathrm{{\mathbb {V}}}\,}}\left[ V_\ell -V_{\ell -1}\right] \) and \({{\,\mathrm{{\mathbb {V}}}\,}}\left[ X_\ell \right] \gg {{\,\mathrm{{\mathbb {V}}}\,}}\left[ V_\ell \right] \), for large \(\Delta t_\ell \).
We first fix \(\Delta t_1 = \varepsilon ^2\) and solve (67) for positive real values of \(M>1\) for which the inequality becomes an equality, taking different values of \(\varepsilon \) and \(t^*\). In Fig. 10 we show these unique, numerically computed M-values. It is clear from the range of the color bar that this M-value consistently exists and lies within the interval [6, 13], meaning that any sign change in the r.h.s. of (67) lies in this range. Next, we evaluate the left hand side of (67) for \(M=6\) and \(M=13\) and plot the results in Fig. 11. From this figure, it is clear that (67) holds for small values of M, but no longer holds once the threshold value from Fig. 10 is passed, meaning that (67) consistently holds once M is sufficiently large. We can thus conclude that it makes sense, from a theoretical point of view, to add an extra coarse level, if it is sufficiently coarse in comparison with the level with \(\Delta t = \varepsilon ^2\).
Rights and permissions
About this article
Cite this article
Løvbak, E., Samaey, G. & Vandewalle, S. A multilevel Monte Carlo method for asymptotic-preserving particle schemes in the diffusive limit. Numer. Math. 148, 141–186 (2021). https://doi.org/10.1007/s00211-021-01201-y
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00211-021-01201-y