A discontinuous Galerkin fast spectral method for multi-species full Boltzmann on streaming multi-processors

TL;DR

This paper introduces a fast spectral discontinuous Galerkin method for solving the multi-species full Boltzmann equation efficiently on GPUs, significantly outperforming traditional stochastic methods in speed.

Contribution

The paper presents a novel GPU-accelerated DGFS implementation for the full Boltzmann equation, achieving high parallel efficiency and speedup over existing stochastic approaches.

Findings

Solves the Boltzmann equation in minutes, two orders faster than DSMC.

Achieves parallel efficiency of 0.96--0.99 on 36 GPUs.

Demonstrates effective weak/strong scaling performance.

Abstract

When the molecules of a gaseous system are far apart, say in microscale gas flows where the surface to volume ratio is high and hence the surface forces dominant, the molecule-surface interactions lead to the formation of a local thermodynamically non-equilibrium region extending few mean free paths from the surface. The dynamics of such systems is accurately described by Boltzmann equation. However, the multi-dimensional nature of Boltzmann equation presents a huge computational challenge. With the recent mathematical developments and the advent of petascale, the dynamics of full Boltzmann equation is now tractable. We present an implementation of the recently introduced multi-species discontinuous Galerkin fast spectral (DGFS) method for solving full Boltzmann on streaming multi-processors. The present implementation solves the inhomogeneous Boltzmann equation in span of few minutes, making it at least two order-of-magnitude faster than the present state-of-art stochastic method---direct simulation Monte Carlo---widely used for solving Boltzmann equation. Various performance metrics, such as weak/strong scaling have been presented. A parallel efficiency of 0.96--0.99 is demonstrated on 36 Nvidia Tesla-P100 GPUs.

Figures

Tables

Table 1. Table 1 . Efficiency and accuracy L ∞ superscript 𝐿 L^{\infty} error ℰ ( p ) = ‖ ∂ t f a ​ n ​ a ​ l ​ y ​ t ​ i ​ c ​ a ​ l ( p ) − ∂ t f n ​ u ​ m ​ e ​ r ​ i ​ c ​ a ​ l ( p ) ‖ , p = { 1 , 2 } formulae-sequence superscript ℰ 𝑝 norm subscript 𝑡 superscript subscript 𝑓 𝑎 𝑛 𝑎 𝑙 𝑦 𝑡 𝑖 𝑐 𝑎 𝑙 𝑝 subscript 𝑡 superscript subscript 𝑓 𝑛 𝑢 𝑚 𝑒 𝑟 𝑖 𝑐 𝑎 𝑙 𝑝 𝑝 1 2 \mathcal{E}^{(p)}=\|\partial_{t}f_{analytical}^{(p)}-\partial_{t}f_{numerical}^{(p)}\|,\;p=\{1,2\} for spatially homogeneous Krook-Wu solution at t = 5.5 𝑡 5.5 t=5.5 for different mass-ratios. N 𝑁 N , N ρ subscript 𝑁 𝜌 N_{\rho} , and M 𝑀 M respectively, denote the number of discretization points in the velocity space, number of Gauss quadrature points in the radial direction, and number of discretization points on full sphere. A fixed velocity domain [ − 12 ,  12 ] 3 superscript 12 12 3 [-12,\;12]^{3} has been used for all the mass-ratios.

$N$	$N_{\rho }$	$m_q/m_p=1$						$m_q/m_p=4$
		$M=6$			$M=12$			$M=6$			$M=12$
		$time\left(s\right)$	${ℰ}^{\left(1\right)}$	${ℰ}^{\left(2\right)}$	$time\left(s\right)$	${ℰ}^{\left(1\right)}$	${ℰ}^{\left(2\right)}$	$time\left(s\right)$	${ℰ}^{\left(1\right)}$	${ℰ}^{\left(2\right)}$	$time\left(s\right)$	${ℰ}^{\left(1\right)}$	${ℰ}^{\left(2\right)}$
16	4	0.00039	3.27e-03	3.27e-03	0.00050	1.77e-03	1.77e-03	0.00039	4.80e-03	1.22e-03	0.00050	3.63e-03	2.47e-04
	8	0.00056	3.73e-03	3.73e-03	0.00085	2.00e-03	2.00e-03	0.00050	4.96e-03	1.33e-03	0.00093	3.62e-03	2.42e-04
	16	0.00084	3.73e-03	3.73e-03	0.00147	2.00e-03	2.00e-03	0.00084	4.96e-03	1.33e-03	0.00148	3.62e-03	2.42e-04
24	6	0.00068	1.37e-04	1.37e-04	0.00114	1.01e-04	1.01e-04	0.00068	1.81e-03	2.06e-02	0.00114	1.79e-03	5.15e-03
	12	0.00117	1.49e-04	1.49e-04	0.00209	9.64e-05	9.64e-05	0.00114	2.12e-03	1.87e-02	0.00210	2.13e-03	6.01e-03
	24	0.00210	1.49e-04	1.49e-04	0.00401	9.64e-05	9.64e-05	0.00210	2.12e-03	1.87e-02	0.00401	2.13e-03	6.01e-03
32	8	0.00159	3.04e-05	3.04e-05	0.00287	2.51e-05	2.51e-05	0.00157	1.54e-04	1.62e-02	0.00286	1.52e-04	1.13e-02
	16	0.00286	3.17e-05	3.17e-05	0.00541	2.45e-05	2.45e-05	0.00286	5.91e-05	1.69e-02	0.00542	5.87e-05	1.03e-02
	32	0.00543	3.17e-05	3.17e-05	0.01057	2.45e-05	2.45e-05	0.00542	5.91e-05	1.69e-02	0.01059	5.87e-05	1.03e-02
40	10	0.00328	1.38e-06	1.38e-06	0.00626	1.26e-06	1.26e-06	0.00326	5.53e-05	4.31e-03	0.00626	5.56e-05	4.35e-03
	20	0.00625	9.35e-07	9.35e-07	0.01226	8.10e-07	8.10e-07	0.00626	5.01e-05	4.29e-03	0.01222	4.97e-05	4.54e-03
	40	0.01227	9.35e-07	9.35e-07	0.02446	8.10e-07	8.10e-07	0.01219	5.01e-05	4.29e-03	0.02431	4.97e-05	4.54e-03
48	12	0.00656	1.04e-07	1.04e-07	0.01289	9.99e-08	9.99e-08	0.00658	8.46e-06	5.76e-04	0.01291	8.45e-06	5.93e-04
	24	0.01290	1.05e-07	1.05e-07	0.02556	9.95e-08	9.95e-08	0.01291	7.17e-06	5.80e-04	0.02561	7.51e-06	6.09e-04
	48	0.02545	1.05e-07	1.05e-07	0.05169	9.95e-08	9.95e-08	0.02550	7.17e-06	5.80e-04	0.05215	7.51e-06	6.09e-04
56	14	0.01204	9.80e-08	9.80e-08	0.02350	9.79e-08	9.79e-08	0.01202	5.22e-06	2.32e-04	0.02352	4.08e-06	1.88e-04
	28	0.02354	9.80e-08	9.80e-08	0.04667	9.79e-08	9.79e-08	0.02353	5.09e-06	2.24e-04	0.04662	3.97e-06	1.87e-04
	56	0.04664	9.80e-08	9.80e-08	0.09303	9.79e-08	9.79e-08	0.04674	5.09e-06	2.24e-04	0.09313	3.97e-06	1.87e-04

Table 2. Table 2 . Numerical parameters for Couette flow (Jaiswal et al . , 2019b ) . Based upon our observations from Table  1 , we have used N ρ = 8 subscript 𝑁 𝜌 8 N_{\rho}=8 , in contrast to N ρ = 32 subscript 𝑁 𝜌 32 N_{\rho}=32 used in (Jaiswal et al . , 2019b ) . This does not affect the recovered bulk properties as illustrated in Fig.  1 , however, it speeds up the computation by a factor of 4.

Parameter	Case C-01
Molecular mass: $\{m_1,m_2\}$ ($\times {10}^{27}kg$)	$\{66.3$, $139.1\}$
Non-dim physical space	$[0,\mathrm{\hspace{0.17em}1}]$
Non-dim velocity space	${[-7,\mathrm{\hspace{0.17em}7}]}^3$
$\{N^3,N_{\rho },M\}$	$\{{32}^3,\mathrm{\hspace{0.17em}8},\mathrm{\hspace{0.17em}12}\}$
Spatial elements	$4$
DG order	3
Time stepping	Euler
Viscosity index: ${\omega }_{\{11,\mathrm{\hspace{0.17em}12},\mathrm{\hspace{0.17em}21},\mathrm{\hspace{0.17em}22}\}}$	$\{0.81,\mathrm{\hspace{0.17em}0.805},\mathrm{\hspace{0.17em}0.805},\mathrm{\hspace{0.17em}0.8}\}$
Scattering parameter: ${\alpha }_{\{11,\mathrm{\hspace{0.17em}12},\mathrm{\hspace{0.17em}21},\mathrm{\hspace{0.17em}22}\}}$	$\{1.4,\mathrm{\hspace{0.17em}1.36},\mathrm{\hspace{0.17em}1.36},\mathrm{\hspace{0.17em}1.32}\}$
Ref. diameter: $d_{\text{ref},pq}$ ($\times {10}^{10}m$)	$\{4.11,\mathrm{\hspace{0.17em}4.405},\mathrm{\hspace{0.17em}4.405},\mathrm{\hspace{0.17em}4.7}\}$
Ref. temperature: $T_{\text{ref},pq}$ ($K$)	$\left\{273\right\}$
Characteristic mass: $m_0$ ($\times {10}^{27}kg$)	$66.3$
Characteristic length: $H_0$ ($mm$)	1
Characteristic velocity: ${𝐮}_0$ ($m/s$)	337.2
Characteristic temperature: $T_0$ ($K$)	273
Characteristic number density: $n_0$ ($m^{-3}$)	$1.680\times {10}^{21}$
Initial conditions
Velocity: $𝐮$ ($m/s$)	0
Temperature: $T$ ($K$)	273
Number density: $n^{\left(1\right)}$ ($m^{-3}$)	$1.680\times {10}^{21}$
Number density: $n^{\left(2\right)}$ ($m^{-3}$)	$8.009\times {10}^{20}$
Knudsen number: $({\mathrm{Kn}}_{11},{\mathrm{Kn}}_{22})$	$(0.793,\mathrm{\hspace{0.17em}0.606})$
Knudsen number: $({\mathrm{Kn}}_{12},{\mathrm{Kn}}_{21})$	$(0.803,\mathrm{\hspace{0.17em}0.555})$
Left wall (purely diffuse) boundary conditions (subscript $l$)
Velocity: ${𝐮}_l$ ($m/s$)	$(0,-50,\mathrm{\hspace{0.17em}0})$
Temperature: $T_l$ ($K$)	273
Right wall (purely diffuse) boundary conditions (subscript $r$)
Velocity: ${𝐮}_r$ ($m/s$)	$(0,+50,\mathrm{\hspace{0.17em}0})$
Temperature: $T_r$ ($K$)	273

Table 3. Table 3 . Performance of the solver for Couette flow test cases. The phase-space is defined using a convenient triplet notation N e /K/N 3 , which corresponds to N e elements in physical space, K order DG (equivalently N p = K − 1 subscript 𝑁 𝑝 𝐾 1 N_{p}=K-1 order polynomial for 1-D domain), and N 3 points in velocity space. n 𝑛 n G ( n > 1 𝑛 1 n>1 ) denotes GPU/CUDA/MPI/parallel execution on n 𝑛 n GPUs shared equally across ( n / 3 ) 𝑛 3 (n/3) nodes. Work units represent the total simulation time for first 52 timesteps. Efficiency is defined as ratio ( 1 1 1 G / n 𝑛 n G )/n, where 1 G and n 𝑛 n G are execution-times on one GPU and n 𝑛 n GPU respectively. M = 12 𝑀 12 M=12 and N ρ = 8 subscript 𝑁 𝜌 8 N_{\rho}=8 is used for all cases.

Phase space	Work Units (s)							Efficiency
	1G	3G	6G	9G	12G	24G	36G	1G/3G	1G/6G	1G/9G	1G/12G	1G/24G	1G/36G
$72/3/{20}^3$	47.580	16.155	8.339	5.698	4.392	2.423	1.774	0.98	0.95	0.93	0.90	0.82	0.84
$72/3/{32}^3$	126.601	42.616	21.551	14.563	11.038	5.784	4.030	0.99	0.98	0.97	0.96	0.91	0.98
$72/3/{48}^3$	391.943	131.081	65.913	44.218	33.513	17.224	11.621	1.00	0.99	0.98	0.97	0.95	1.05
$72/6/{20}^3$	94.682	31.957	16.197	10.944	8.331	4.392	3.079	0.99	0.97	0.96	0.95	0.90	0.96
$72/6/{32}^3$	253.016	84.834	42.741	28.697	21.703	11.158	7.693	0.99	0.99	0.98	0.97	0.94	1.03
$72/6/{48}^3$	782.343	261.601	131.217	87.755	66.009	33.520	22.509	1.00	0.99	0.99	0.99	0.97	1.09
$216/3/{20}^3$	141.754	47.641	24.033	16.182	12.326	6.356	4.388	0.99	0.98	0.97	0.96	0.93	1.01
$216/3/{32}^3$	378.956	126.853	63.676	42.636	32.066	16.295	11.041	1.00	0.99	0.99	0.98	0.97	1.07
$216/3/{48}^3$	1172.907	391.916	196.439	131.153	98.538	49.652	33.471	1.00	1.00	0.99	0.99	0.98	1.10
$216/6/{20}^3$	283.091	94.737	47.679	31.903	24.060	12.262	8.320	1.00	0.99	0.99	0.98	0.96	1.06
$216/6/{32}^3$	759.149	253.498	127.004	84.932	63.780	32.212	21.672	1.00	1.00	0.99	0.99	0.98	1.09
$216/6/{48}^3$	2347.099	783.642	392.470	261.817	196.552	98.680	66.018	1.00	1.00	1.00	1.00	0.99	1.11

Equations

$$\frac{\partial f^{(p)}}{\partial t}+\mathbf{c}\cdot {\nabla }_{\hat{\mathbf{x}}}{f}^{(p)}=\sum _q\frac{1}{{\mathrm{Kn}}_{pq}}{\mathcal{Q}}^{(pq)},p=1,2,\dots ,n,$$

$${\mathcal{Q}}^{(pq)}(f^{(p)},f^{(q)})(\mathbf{c})$$

$$\left[f^{(p)}({\mathbf{c}}^{\prime })f^{(q)}({\mathbf{c}}_{\ast }^{\prime })-f^{(p)}(\mathbf{c})f^{(q)}({\mathbf{c}}_{\ast })\right]\mathrm{d}\mathbf{\sigma }\mathrm{d}{\mathbf{c}}_{\ast },$$

$$\displaystyle\left\{\begin{array}[]{l}\displaystyle\mathbf{c}^{\prime}=\frac{\mathbf{c}+\mathbf{c}_{*}}{2}+\frac{(1-m_{q}/m_{p})}{2(1+m_{q}/m_{p})}(\mathbf{c}-\mathbf{c}_{*})+\frac{1}{(m_{p}/m_{q}+1)}|\mathbf{c}-\mathbf{c}_{*}|{\boldsymbol{\sigma}},\\[12.0pt] \displaystyle\mathbf{c}_{*}^{\prime}=\frac{\mathbf{c}+\mathbf{c}_{*}}{2}+\frac{(1-m_{q}/m_{p})}{2(1+m_{q}/m_{p})}(\mathbf{c}-\mathbf{c}_{*})-\frac{1}{(m_{q}/m_{p}+1)}|\mathbf{c}-\mathbf{c}_{*}|{\boldsymbol{\sigma}},\end{array}\right.$$

$$B_{(pq)}$$

$$\frac{{\alpha }_{pq}}{2^{1+{\alpha }_{pq}}\Gamma (2.5-{\omega }_{pq})\pi }|\mathbf{c}-{\mathbf{c}}_{\ast }{|}^{2(1-{\omega }_{pq})}(1+\cos \chi {)}^{{\alpha }_{pq}-1},$$

$${\mathrm{Kn}}_{pq}$$

$${\beta }^{(pq)}$$

$${\mathcal{Q}}^{(pq)}(f^{(p)},f^{(q)})(\mathbf{c})$$

$$\displaystyle\Bigl{[}f^{(p)}(\mathbf{c}^{\prime})f^{(q)}(\mathbf{c}^{\prime}_{*})-f^{(p)}(\mathbf{c})f^{(q)}(\mathbf{c}-\mathbf{u})\Bigr{]}\mathrm{d}{{\boldsymbol{\sigma}}}\mathrm{d}{\mathbf{u}},$$

$$\displaystyle\left\{\begin{array}[]{l}\displaystyle\mathbf{c}^{\prime}=\mathbf{c}-\frac{m_{q}}{m_{p}+m_{q}}\mathbf{u}+\frac{m_{q}}{m_{p}+m_{q}}|\mathbf{u}|{\boldsymbol{\sigma}},\\[12.0pt] \displaystyle\mathbf{c}_{*}^{\prime}=\mathbf{c}-\frac{m_{q}}{m_{p}+m_{q}}\mathbf{u}-\frac{m_{p}}{m_{p}+m_{q}}|\mathbf{u}|{\boldsymbol{\sigma}}.\end{array}\right.$$

$$R=\frac{4}{1+\max (4m_q/(m_p+m_q),2)+\sqrt{1+m_q/m_p}}L$$

$$f^{(p)}(\mathbf{c})=\sum _{k=-N/2}^{N/2-1}{\hat{f}}_k^{(p)}{e}^{i\frac{\pi }{L}k\cdot \mathbf{c}},f^{(q)}(\mathbf{c})=\sum _{k=-N/2}^{N/2-1}{\hat{f}}_k^{(q)}{e}^{i\frac{\pi }{L}k\cdot \mathbf{c}}.$$

$$\begin{split}\hat{\mathcal{Q}}^{(pq)}_{k}:&=\frac{1}{(2L)^{3}}\int_{D_{L}}\mathcal{Q}^{(pq)}(f^{(p)},f^{(q)})(\mathbf{c})e^{-i\frac{\pi}{L}k\cdot\mathbf{c}}\mathrm{d}{\mathbf{c}}\\ &=\sum_{\begin{subarray}{c}l,m=-N/2\\ l+m=k\end{subarray}}^{N/2-1}\Bigl{[}G^{{(pq)}+}(l,m)-G^{{(pq)}-}(m,m)\Bigr{]}\;\hat{f}^{(p)}_{l}\hat{f}^{(q)}_{m},\end{split}$$

$$G^{(pq)+}(l,m)$$

$$\left[e^{-i\frac{\pi }{L}\frac{m_q}{m_p+m_q}(l+m)\cdot \mathbf{u}+i\frac{\pi }{L}|\mathbf{u}|\left(\frac{m_q}{m_p+m_q}l-\frac{m_p}{m_p+m_q}m\right)\cdot \mathbf{\sigma }}\right]\mathrm{d}\mathbf{\sigma }\mathrm{d}\mathbf{u}$$

$$G^{(pq)-}(m,m)$$

$$G^{(pq)+}(l,m)\approx \sum _{r=1}^{N_{\rho }}{\alpha }_r(l+m){\beta }_r(l){\gamma }_r(m),$$

$${\hat{\mathcal{Q}}}_k^{(pq)+}=\sum _{r=1}^{N_{\rho }}{\alpha }_r(k)\sum _{\begin{array}{c}l,m=-N/2\\ l+m=k\end{array}}^{N/2-1}\left({\beta }_r(l){\hat{f}}_l^{(p)}\right)\left({\gamma }_r(m){\hat{f}}_m^{(q)}\right),$$

$$G^{(pq)+}(l,m)=$$

$$e^{-i\frac{\pi }{L}\frac{m_q}{m_p+m_q}(l+m)\cdot \mathbf{u}+i\frac{\pi }{L}|\mathbf{u}|\left(\frac{m_q}{m_p+m_q}l-\frac{m_p}{m_p+m_q}m\right)\cdot \mathbf{\sigma }}\mathrm{d}\mathbf{\sigma }\mathrm{d}\mathbf{u}$$

$$={\int }_0^R{\int }_{S^{d-1}}$$

$$F^{(pq)}(l+m,\rho ,\mathbf{\sigma })={\rho }^2{\int }_{S^2}{B}_{pq}(\rho ,\mathbf{\sigma }\cdot \hat{\mathbf{u}})e^{-i\frac{\pi }{L}\rho \frac{m_q}{m_p+m_q}(l+m)\cdot \hat{\mathbf{u}}}\mathrm{d}\hat{\mathbf{u}},$$

$$G^{(pq)-}(m)$$

$$={\int }_0^R{\int }_{S^2}{\int }_{S^2}{\rho }^2{B}_{pq}(\rho ,\mathbf{\sigma }\cdot \hat{\mathbf{u}})e^{-i\frac{\pi }{L}\rho m\cdot \hat{\mathbf{u}}}\mathrm{d}\mathbf{\sigma }\mathrm{d}\hat{\mathbf{u}}\mathrm{d}\rho .$$

$${\partial }_t{f}^{(p)}=\sum _{q=1}^2{\int }_{{\mathbb{R}}^3}{\int }_{S^2}{B}_{pq}\left[f^{(p)}(v^{\prime })f^{(q)}(v_{\ast }^{\prime })-f^{(p)}(v)f^{(q)}(v_{\ast })\right]\mathrm{d}\mathbf{\sigma }\mathrm{d}{v}_{\ast },$$

$$\displaystyle f^{(p)}(t,v)=n^{(p)}\Bigg{(}\frac{m_{p}}{2\pi K}\Bigg{)}^{3/2}\exp{\Bigg{(}-\frac{m_{p}v^{2}}{2K}\Bigg{)}}\Bigg{(}(1-3Q_{p})+\frac{m_{p}}{K}Q_{p}v^{2}\Bigg{)},$$

$$\mu =\frac{4m_1{m}_2}{(m_1+m_2{)}^2},{\tau }_1={\lambda }_{22}-{\lambda }_{21}\mu (3-2\mu ),{\tau }_2={\lambda }_{11}-{\lambda }_{12}\mu (3-2\mu ),$$

$$\displaystyle A=\frac{1}{6}\Bigg{(}\lambda_{11}+\lambda_{21}\mu\left(3-2\mu\frac{\tau_{2}}{\tau_{1}}\right)\Bigg{)},\quad B=\frac{1}{3}\Bigg{(}\lambda_{11}\tau_{1}+\lambda_{21}\mu(3-2\mu)\tau_{2}\Bigg{)},$$

$$Q(t)=\frac{A}{A\exp (At)-B},Q_p(t)={\tau }_{p}Q(t),$$

$$K(t)=\frac{n^{(1)}+n^{(2)}}{(n^{(1)}+n^{(2)})+2(n^{(1)}{\tau }_1+n^{(2)}{\tau }_2)Q(t)}.$$

$$({\tau }_1-{\tau }_2)\left(2{\mu }^2\left(\frac{{\lambda }_{21}}{{\tau }_1}-\frac{{\lambda }_{12}}{{\tau }_2}\right)-1\right)=0.$$

$$\displaystyle\partial_{t}f^{(p)}=f^{(p)}\Bigg{(}-\frac{3}{2K}K^{\prime}+\frac{m_{p}\;v^{2}}{2K^{2}}K^{\prime}\Bigg{)}$$

$$\displaystyle+n^{(p)}\Bigg{(}\frac{m_{p}}{2\pi K}\Bigg{)}^{3/2}\exp{\Bigg{(}-\frac{m_{p}v^{2}}{2K}\Bigg{)}}\Bigg{(}-3Q^{\prime}_{p}+\frac{m_{p}}{K}Q^{\prime}_{p}v^{2}-\frac{m_{p}}{K^{2}}K^{\prime}Q_{p}v^{2}\Bigg{)}$$

$$:=\sum _{q=1}^2{\mathcal{Q}}^{(pq)}(f^p,f^q),$$

$$Q^{\prime }(t)=-\frac{A^3\exp (At)}{(A\exp (At)-B{)}^2},Q_p^{\prime }(t)={\tau }_p{Q}^{\prime }(t),$$

$$K^{\prime }(t)=-\frac{2(n^{(1)}+n^{(2)})(n^{(1)}{\tau }_1+n^{(2)}{\tau }_2)}{[(n^{(1)}+n^{(2)})+2(n^{(1)}{\tau }_1+n^{(2)}{\tau }_2)Q(t){]}^2}{Q}^{\prime }(t).$$

Full text

Multi-species full Boltzmann, Discontinuous Galerkin Fast Spectral, Graphics Processing Units

††copyright: rightsretained††conference: Platform for Advanced Scientific Computing; June 2019; ETH Zurich, Switzerland††ccs: Applied computing Mathematics and statistics††ccs: Computing methodologies Massively parallel and high-performance simulations

1. Introduction

From the fundamental mass/momentum conservation principles, it can be inferred that, in the presence of external forces, say, pressure and temperature gradients, the heavier species moves slower and the lighter species moves faster giving rise to a phenomena termed as diffusion, and effects thereof. Diffusion processes are critical in many applications, for instance, the measurement of the neutrino mass using a windowless gaseous tritium source in the ongoing KATRIN experiment (Bornschein et al., 2005). The dynamics of such systems (and others) are governed by the Boltzmann equation—an integro-differential equation describing the evolution of the distribution function in six-dimensional phase space—which models the dilute gas behavior at the molecular level to accurately describe a wide range of non-continuum flow phenomena, for instance, shocks, expansions into vacuum (Muntz, 1989) as well as velocity and thermal slip at gas-solid interfaces (Sharipov and Kalempa, 2003, 2004). Most rarefied flows of technological interest involve gas mixtures with species diffusion playing a decisive role in turbulent, chemically reacting flows, and evaporation/condensation processes (Takata and Golse, 2007).

The approaches for numerical solution of the Boltzmann equation date back to as early as 1940s (Grad, 1949) using, for example, the now widely used direct simulation Monte Carlo (DSMC) method (Bird, 1994). The DSMC method, based on the kinetic theory of dilute gases, models the binary interactions between particles stochastically. However, it is this stochastic nature, that makes the method unsuitable for flows involving species in trace concentration, for instance, to analyze the spectrum of beta electrons emitted by tritium source which can be substantially different in the presence of the impurities in KATRIN experiment (Bornschein et al., 2005; Sharipov and Kalempa, 2005).

The main difficulty of numerically solving the full Boltzmann equation lies in its complicated collision term. Recently, a fast Fourier spectral method for the multi-species Boltzmann collision operator was introduced in (Jaiswal et al., 2019b). The complexity for a single evaluation of the collision operator is reduced from $O(N^{6})$ (direct calculation) to $O(MN_{\rho}N^{3}\log N)$ (based on a low-rank decomposition strategy), where $N$ is the number of discretization points in each velocity dimension, $N_{\rho}\sim O(N)$ is the number of discretization points in the radial direction needed for low-rank decomposition, and $M\ll N^{2}$ is the number of discretization points on the sphere. Based on (Jaiswal et al., 2019a), a discontinuous Galerkin fast spectral (DGFS) method was also proposed in (Jaiswal et al., 2019b) for solving the full multi-species Boltzmann equation. DGFS can produce high order spatially and temporally accurate solutions for low-speed and unsteady flows in micro-systems, and is amenable to excellent nearly-linear scaling characteristics on massively parallel architectures. This paper focuses on implementation aspects of multi-species DGFS with an emphasis on establishing the algorithmic behavior of such numerical schemes. More specifically, here, we are concerned about the scaling characteristics of DGFS on multi-GPU/multi-CPU systems.

In the section that follows, we describe the multi-species Boltzmann in brief, followed by a description of the collision operator algorithm. Various performance metrics such as weak/strong scaling, and micro benchmarks involving, static adaptivity of cell-size and polynomial order approximation for rarefied gas-flows have been discussed in section 3. Concluding remarks are given in section 4.

2. The multi-species Boltzmann equation

The non-dimensional Boltzmann equation for multi-species, mono-atomic gas without external forces can be written as (cf. (Jaiswal et al., 2019b))

[TABLE]

where $n$ denotes number of species in the mixture – each of them represented by a number distribution function $f^{(p)}(t,\mathbf{x},\mathbf{c})$ of time $t$, position $\mathbf{x}$, and particle velocity $\mathbf{c}$. The collision operator $\mathcal{Q}^{(pq)}$ takes into account interactions between species $p$ and $q$, which acts only in the velocity space:

[TABLE]

where $(\mathbf{c},\mathbf{c}_{*})$ and $(\mathbf{c}^{\prime},\mathbf{c}^{\prime}_{*})$ denote the pre and post collision velocity pairs, which are related through momentum and energy conservation as

[TABLE]

where $m_{p}$, $m_{q}$ denote the mass of particles of species $p$ and $q$ respectively. Here, the vector ${\boldsymbol{\sigma}}$ varies over the unit sphere $\mathcal{S}^{2}$. The quantity $\mathcal{B}_{pq}$ ($\geq 0$) is the collision kernel depending only on $|\mathbf{c}-\mathbf{c}_{*}|$ and the scattering angle $\chi$ (angle between $\mathbf{c}-\mathbf{c}_{*}$ and $\mathbf{c}^{\prime}-\mathbf{c}^{\prime}_{*}$). In the present work, we consider the variable soft sphere (VSS) (Koura and Matsumoto, 1991) scattering model. It is worth emphasizing that although the VSS collision kernel is adopted in the present work for easy comparison with DSMC solutions, the fast spectral method we use for the collision operator applies straightforwardly to general collision kernels (see (Gamba et al., 2017; Jaiswal et al., 2019a, b)).

For Variable Soft-Sphere model (Bird, 1994) in particular, the non-dimensional collision-kernel $B_{(pq)}$, and the Knudsen number $\mathrm{Kn}_{pq}$ are given as

[TABLE]

[TABLE]

Here $\Gamma$ denotes the usual Gamma function, $d_{(\mathrm{ref},{pq})}$, $T_{(\mathrm{ref},{pq})}$, $\omega_{pq}$, and $\alpha_{pq}$ are, respectively, the reference diameter, the reference temperature, the viscosity index, and the scattering parameter. The diameter $d_{(\mathrm{ref},{pq})}$ and exponent $\alpha_{pq}$ are determined so that the transport (viscosity and diffusion) coefficients of VSS are consistent with experimental data. Additionally $H_{0}$, $T_{0}$, $n_{0}$, and $m_{0}$, respectively, denote the characteristic length, characteristic temperature, characteristic number density, and characteristic mass $m_{0}$. Based upon these, we define the characteristic velocity as $u_{0}=\sqrt{2k_{B}T_{0}/m_{0}}$ where $k_{B}$ refers to Boltzmann constant; and characteristic time as $t_{0}=H_{0}/u_{0}$. For convenience, we define a pre-factor $\beta^{(pq)}$ as

[TABLE]

Henceforth, we will always refer to the non-dimensional Boltzmann equation (1) in our presentation.

2.1. The collision operator

First, note that $\mathcal{Q}^{(pq)}(f^{(p)},f^{(q)})$ does not depend on spatial coordinate $\mathbf{x}$. Given distribution functions $f^{(p)}$ and $f^{(q)}$ of species $p$ and $q$, dependent only on the velocity coordinate $\mathbf{c}$: discretized uniformly using $N^{3}$ points, the method produces $\mathcal{Q}^{(pq)}(f^{(p)},f^{(q)})$ at the same grid with $O(MN_{\rho}N^{3}\log N)$ complexity, where $N_{\rho}\sim O(N)$ is the number of Gauss-Legendre quadrature/discretization points in the radial direction needed for low-rank decomposition , $M\ll N^{2}$ is the number of discretization points on the sphere. The steps (based on (Jaiswal et al., 2019b)) for evaluating $\mathcal{Q}^{(pq)}$ can be summarized as:

• •

Change the variable $\mathbf{c}_{*}$ to $\mathbf{u}=\mathbf{c}-\mathbf{c}_{*}$:

[TABLE]

where $\hat{\mathbf{u}}$ is the unit vector along $\mathbf{u}$, and

[TABLE]

• •

Determine the extent of velocity domain $D_{L}=[-L,L]^{3}$, and periodically extend $f$, $g$ to $\mathbb{R}^{3}$.

• •

Truncate the integral in $\mathbf{u}$ to a ball $B_{R}$ with

[TABLE]

• •

Approximate $f^{(p)}$, $f^{(q)}$ by truncated Fourier series

[TABLE]

Note here $k$ is a three-dimensional index.

• •

Substitute $f^{(p)}$, $f^{(q)}$ into (9), and perform the standard Galerkin projection

[TABLE]

where $k=-N/2,\dots,N/2-1$, and the kernel modes $G^{{(pq)}+}$ and $G^{{(pq)}-}$ are given by

[TABLE]

It is clear that the direct evaluation of $\hat{\mathcal{Q}}^{(pq)}_{k}$ (for all $k$) would require $O(N^{6})$ complexity. But if we can find a low-rank, separated expansion of $G^{{(pq)}+}(l,m)$ as

[TABLE]

then the gain term (positive part) of $\hat{\mathcal{Q}}^{(pq)}_{k}$ can be rearranged as

[TABLE]

which is a convolution of two functions $\beta_{r}(l)\hat{f}^{(p)}_{l}$ and $\gamma_{r}(m)\hat{f}^{(q)}_{m}$, hence can be computed via fast Fourier transform (FFT) in $O(N_{\rho}N^{3}\log N)$ operations. Note that the loss term (negative part) of $\hat{\mathcal{Q}}^{{(pq)}-}_{k}$ is readily a convolution and can be computed via FFT in $O(N^{3}\log N)$ operations.

In order to find the approximation in (17), we simplify (16) as

[TABLE]

where

[TABLE]

while for the loss term,

[TABLE]

For details on the error introduced from the Fourier-spectral approximation, the reader is referred to (Jaiswal et al., 2019b). This discussion has been omitted in the present work for brevity.

2.2. The collision operator algorithm

The collision operator procedure described above is applicable for general collision kernels for $n$-species mixture. However, for a concise description of the algorithmic ideas from an implementation viewpoint, we restrict our discussion to Variable Soft Sphere collision kernel (2). The ideas, however, can certainly be carried over to other collision kernels.

In multi-species implementation, with the high amount of involved computation, our motive is to avoid spurious computation for every timestep. We first outline the procedure for pre-computing variables that can be stored and reused during the course of the simulation.

• •

First, we precompute $(\pi/L\;\rho\;l\cdot{\boldsymbol{\sigma}})$. We use Gauss-Legendre-Quadrature (GLQ) for integration. So $\rho$, the GLQ zeros, is an array of size $N_{\rho}$ (since the integrand oscillates on the scale of O(N), the total number of quadrature points needed should be $\sim O(N)$). Additionally, we use spherical design (Womersley, 2016) quadrature on sphere. So, ${\boldsymbol{\sigma}}$, the spherical-quadrature zeros, is an array of size $M$. $l$ as previously defined is the 3-D velocity-space index, and is therefore an array of size $N^{3}$. Based upon these $(\pi/L\;\rho\;l\cdot{\boldsymbol{\sigma}})$ is precomputed and stored as a $N_{\rho}\times M\times N^{3}$ flattened row-major array $\mathrm{a}_{xyz}$. This is described in steps 1–9 of Algo. (1).

• •

Second, we compute $F(l+m,\rho,{\boldsymbol{\sigma}})$ as per Eq. (20). Note that $k=l+m$ is velocity-space index of size $N^{3}$. Since $l+m$, $\rho$, and ${\boldsymbol{\sigma}}$ do not change with time, the term $F(l+m,\rho,{\boldsymbol{\sigma}})$ is precomputed and stored as a $N_{\rho}\times M\times N^{3}$ flattened row-major array $\mathrm{b}^{(pq)}_{xyz}$ for every collision pair $(p,q)$. This is described in step 13 of Algo. (1).

• •

Third, we perform precomputation needed for loss-term $G^{{(pq)}-}(m)$ as per Eq. (21). The output is stored as a $N^{3}$ flattened row-major array $\mathrm{c}^{(pq)}_{z}$ for every collision pair $(p,q)$. This is described in step 14 of Algo. (1).

Next we outline the procedure for computing $\mathcal{Q}^{(pq)}$. Recall that our motive is to compute (15)

• •

First, we compute the forward Fourier transform of $\mathcal{F}^{(p)}_{i,\,l_{1}}$, and $\mathcal{F}^{(q)}_{i,\,l_{2}}$ to obtain $\hat{f}^{(p)}_{l}$ and $\hat{f}^{(q)}_{m}$ respectively. This is described in step 1 of Algo. (2).

• •

Second, we compute $G^{{(pq)}+}(l,m)$ as per Eq. (2.1). Recall that $(\pi/L\;\rho\;l\cdot{\boldsymbol{\sigma}})$ has been already precomputed and stored as $\mathrm{a}_{xyz}$. Also recall that $F(l+m,\rho,{\boldsymbol{\sigma}})$ has been precomputed and stored as $\mathrm{b}^{(pq)}_{xyz}$. These can be reused to compute $G(l,m)$. This is described in step 2–8 of Algo. (2). In our implementation, we explicitly unroll the nested loops using Mako (Bayer, 2018) templating engine, such that variables $\mathtt{t1,t2}$ in steps 4 and 5 are computed in a single kernel call (thereby requiring a space of $MN_{\rho}N^{3}$ each), and the FFT transforms in the step 6 are rather $MN_{\rho}$ batched FFT transforms, each of size $N^{3}$.

• •

Third, in order to perform convolution for the loss-term $G^{{(pq)}-}(l,m)$, we prepare the variable QG in step 7 of Algo. (2).

• •

Fourth, we perform convolutions to compute $\hat{\mathcal{Q}}^{(pq)}_{k}$ as in Eq. (15). Recall that $G^{{(pq)}-}(m)$ has now been precomputed and stored as QG, and can be reused here. An inverse Fourier transform is then performed to obtain final $\mathcal{Q}^{(pq)}$. This is described in step 10 of Algo. (2).

3. Micro-Benchmarks

Verification for standard rarefied gas flows can be found in (Jaiswal et al., 2019b). In the present work, we focus on the evaluation of the algorithmic behavior.

3.1. Hardware Configuration

Serial and parallel implementations of multi-species DGFS solver are run on 15-node Brown-GPU RCAC cluster at Purdue University. Each node is equipped with two 12-core Intel Xeon Gold 6126 CPU, and three Tesla-P100 GPU. The operating system used is 64-bit CentOS 7.4.1708 (Core) with NVIDIA Tesla-P100 GPU accompanying CUDA driver 8.0 and CUDA runtime 8.0. The GPU has 10752 CUDA cores, 16GB device memory, and compute capability of 6.0. The solver has been written in Python/PyCUDA and is compiled using OpenMPI 2.1.0, g++ 5.2.0, and nvcc 8.0.61 compiler with third level optimization flag. All the simulations are done with double precision floating point values.

3.2. Spatially homogeneous case: Krook-Wu exact solution

For constant collision kernel, an exact solution to the spatially homogeneous multi-species Boltzmann equation can be constructed (see (Krook and Wu, 1977)). We use this solution to verify the accuracy of the proposed fast spectral method for approximating the collision operator. Considering a binary mixture ($n=2$; $p=1,2$), the equation simplifies to

[TABLE]

where $B_{{pq}}=B_{{qp}}:=\frac{\lambda_{{qp}}}{4\pi n^{{(q)}}}$ and $\lambda_{{pq}}$ is some positive constant. The exact solution is given by

[TABLE]

where

[TABLE]

Furthermore, the following condition needs to be satisfied

[TABLE]

For simplicity, we choose $n^{(1)}=n^{(2)}=1$, $\lambda_{11}=\lambda_{22}=1$, $\lambda_{12}=\lambda_{21}=1/2$ but vary the mass ratio $m_{1}/m_{2}$ in the following tests.

It is also helpful to take the derivative of eqn. (23), which yields

[TABLE]

where

[TABLE]

This allows us to check the accuracy of the collision solver without introducing time discretization error.

Table 1 shows the $L^{\infty}$ norm between the numerical and analytical $\partial f^{(p)}/\partial t$. For different mass ratios, we have considered the cases with $N=\{16,\,24,\,32,\,40,\,48,\,56\}$ points in each velocity dimension; and $M=6$, $12$ spherical design quadrature points on the full sphere. A good agreement between analytical and numerical solutions is evident from the table. At a fixed $N$, with increase in mass ratio, the error norm increases. In particular, increase in $M$ does not considerably affect the solution due to the isotropic nature of the distribution function. Note that, in the fast spectral decomposition, since the integral oscillates roughly on $O(N)$, the total number of Gauss–Legendre quadrature points $N_{\rho}$ in the radial direction should be on order of $O(N)$. As per (Gamba et al., 2017), a more precise estimate is $\approx 0.8N$. However, there is no good rule to select optimal $N_{\rho}$. We observe that the error is relatively unaffected upon reducing $N_{\rho}$ from $N$ to $N/2$. However, we note that $N_{\rho}=N$ is a safer choice.

From a computational viewpoint, the simulation time is independent of the mass ratio. On increasing the number of discretization points on the sphere $M$, the computational cost approximately doubles–however, we do observe the effect of loop unrolling for smaller $N$. Likewise, the computational cost approximately doubles on increasing the number of quadrature points $N_{\rho}$. This establishes that the algorithm is linear in both $M$ and $N_{\rho}$.

3.3. Spatially in-homogeneous case: Couette flow

The aforementioned methodology allows us to compute the collision operator efficiently. To solve the fully spatial in-homogeneous equation (1), we also need an accurate and efficient spatial and time discretization. Here, we adopt the the Runge-Kutta discontinuous Galerkin (RKDG) approach–widely used for hyperbolic systems–as adapted in (Jaiswal et al., 2019a, b) for Boltzmann equation. The details of the discretization can be found in (Jaiswal et al., 2019a, b). We mention that evaluation of collision operator consumes $>98\%$ of computation time, and hence, in the present work, we focus on the collision operator behavior. More details on spatial-temporal RKDG discretization on GPU can be found in (Klöckner et al., 2009; Witherden et al., 2014). We restrict our discussion and benchmarks to 1-D flow problems for brevity111Discussion and benchmarks for higher 2D/3D spatial dimension shall be presented in the extended version of this manuscript..

3.3.1. Verification

For general Boltzmann equation (1), analytical solutions do not exist. Therefore, we compare our results with widely accepted direct simulation Monte Carlo (DSMC) (Bird, 1994) method. We want to emphasize that DSMC is a stochastic method for solution of the N-particle master kinetic equation which converges to the Boltzmann equation in the limit of infinite number of particles (Wagner, 1992).

In the current test case, we consider the effect of velocity gradient on the solution. The coordinates are chosen such that the walls are parallel to the $y$ direction and $x$ is the direction perpendicular to the walls. The geometry as well as boundary conditions are shown in Figure 1. Figure 2 illustrates the velocity and temperature along the domain length for both species, wherein we observe an excellent agreement between DGFS and DSMC. The small discrepancies, however, are primarily due to: a) statistical fluctuations inherent to the Monte Carlo methods, b) practical limitations on number of particles used in DSMC simulations. From a computational viewpoint, the present DGFS simulations on a single GPU took 138 seconds to acquire the steady state, in contrast to 26086.45 sec on 24 processors for DSMC simulations as reported in (Jaiswal et al., 2019b), for achieving comparable accuracy.

3.3.2. Scaling Behavior

The simulations are carried out for different test-cases by varying element-count ($N_{e}$), polynomial approximation order ($N_{p}=K-1$), and velocity-space sizes ($N$). The spatial elements are distributed to $p$ processors using the well-known linear domain-decomposition strategy requiring sharing of $O(pN^{3})$ floating-point during MPI communication phase. Speed up obtained with multi-GPU solver is presented in Table (3). As evident from the table, the acceleration due to GPU parallelization increases with increase in the size of computational grid. More specifically, the increase in $N_{e}$ and $K$ have small-effect on overall speedup which suggests that DG-operators (for instance derivative, time-evolution) are rather computationally inexpensive operations. On the other hand, increase in velocity-grid improves the observed speedup. The weak/strong scaling behavior is also evident from the table.

3.3.3. Flat profile

Recall that the fast Fourier spectral collision operator algorithm 2 is split into multiple parts. It is therefore interesting to see what performance level is attained by each part of the operator. Fig (3) presents the percentage of time spent in various parts of Algo. 2 vs. order of DG scheme (K). First, we note that the DG operators denoted in yellow, requires $~{}1\%$ of the total simulation time. The collision operator, however, consumes nearly $>98\%$ of the total time for both $N^{3}=20^{3}$ and $N^{3}=32^{3}$.

4. Conclusions

We have presented an implementation of the multi-species Discontinuous Galerkin Fast Spectral (DGFS) method for solution of multi-species monoatomic full Boltzmann equation on multi-GPU/multi-CPU architectures. The DG-type formulation employed in the present work has advantage of having high-order accuracy at the element-level, and its element-local compact nature (and that of our collision algorithm) enables effective parallelization on massively parallel architectures. For verification and benchmarks, we carry out simulations for spatially homogeneous BKW, and Couette flow problems. Parallel efficiency close to 0.95 is observed on a 36 GPU multi-node/multi-GPU system. An important key observation is that the efficiency can be maintained provided we have enough work on each processor. It is this speedup that now allows researchers to solve problems within a day that would otherwise take months on traditional CPUs. Future work directions include, assessment of the implementation beyond thousand cores. Extending the implementation to general 2D/3D mixed grids coupled with adaptivity in physical and velocity spaces, is an interesting direction as well.