Fast Kinetic Scheme : efficient MPI parallelization strategy for 3D Boltzmann equation

TL;DR

This paper introduces an MPI parallelization strategy for the Fast Kinetic Scheme to efficiently solve the 3D Boltzmann equation, demonstrating near-ideal scalability on large-scale systems.

Contribution

It extends the Fast Kinetic Scheme to the full Boltzmann equation and develops an efficient parallelization approach for high-dimensional kinetic simulations.

Findings

Achieves close to ideal scalability on tera- and peta-scale systems.

Successfully extends the scheme from BGK approximation to full Boltzmann equation.

Demonstrates efficient parallelization for high-dimensional kinetic equations.

Abstract

In this paper we present a parallelization strategy on distributed memory systems for the Fast Kinetic Scheme --- a semi-Lagrangian scheme developed in [J. Comput. Phys., Vol. 255, 2013, pp 680-698] for solving kinetic equations. The original algorithm was proposed for the BGK approximation of the collision kernel. In this work we deal with its extension to the full Boltzmann equation in six dimensions, where the collision operator is resolved by means of fast spectral method. We present close to ideal scalability of the proposed algorithm on tera- and peta-scale systems.

Figures

Tables

Table 1. Table 1: Computational time, number of neighboring nodes, number of data cells per node and number of ghost cells per node as a function of MPI size ( N M ​ P ​ I subscript 𝑁 𝑀 𝑃 𝐼 N_{M\!P\!I} ) and domain decomposition (MPI dimensions N M ​ P ​ I x subscript 𝑁 𝑀 𝑃 subscript 𝐼 𝑥 N_{M\!P\!I_{x}} , N M ​ P ​ I y subscript 𝑁 𝑀 𝑃 subscript 𝐼 𝑦 N_{M\!P\!I_{y}} , N M ​ P ​ I z subscript 𝑁 𝑀 𝑃 subscript 𝐼 𝑧 N_{M\!P\!I_{z}} ) for the BGK kernel test case on the 64 3 × 32 3 superscript 64 3 superscript 32 3 64^{3}\times 32^{3} mesh. The smallest number of ghost cells and the smallest run time are written in bold type for every tested MPI size N M ​ P ​ I subscript 𝑁 𝑀 𝑃 𝐼 N_{M\!P\!I} .

$N_{MPI}$	$N_{MP{I}_x}$	$N_{MP{I}_y}$	$N_{MP{I}_z}$	time	#neigh.	#cells	#ghosts
2	2	1	1	$1376.51$	2	131072	8192
4	4	1	1	$704.064$	2	65536	8192
	2	2	1	$712.356$	8	65536	8448
8	8	1	1	$359.652$	2	32768	8192
	4	2	1	$362.894$	8	32768	6400
	2	2	2	$386.392$	26	32768	6536
16	16	1	1	$191.856$	2	16384	8192
	8	2	1	$189.828$	8	16384	5376
	4	4	1	$191.113$	8	16384	4352
	4	2	2	$201.823$	26	16384	4424
32	32	1	1	$106.414$	2	8192	8192
	16	2	1	$101.946$	8	8192	4864
	8	4	1	$100.669$	8	8192	3328
	8	2	2	$106.942$	26	8192	3368
	4	4	2	$108.903$	26	8192	2824
64	64	1	1	$63.5180$	2	4096	8192
	32	2	1	$58.7658$	8	4096	4608
	16	4	1	$55.6507$	8	4096	2816
	16	2	2	$59.4585$	26	4096	2840
	8	8	1	$56.5515$	8	4096	2304
	8	4	2	$59.3814$	26	4096	2024
	4	4	4	$62.6682$	26	4096	1736

Table 2. Table 2: Performance tests on 270 3 × 16 3 superscript 270 3 superscript 16 3 270^{3}\times 16^{3} mesh for BGK. Time per cycle is obtained by T cycle = T / N cycle subscript 𝑇 cycle 𝑇 subscript 𝑁 cycle T_{\text{cycle}}=T/N_{\text{cycle}} , time per cycle per cell by T cell = T cycle / N c subscript 𝑇 cell subscript 𝑇 cycle subscript 𝑁 𝑐 T_{\text{cell}}=T_{\text{cycle}}/N_{c} and time per cycle per node by T cycle/node = N s ​ T cell subscript 𝑇 cycle/node subscript 𝑁 𝑠 subscript 𝑇 cell T_{\text{cycle/node}}=N_{s}T_{\text{cell}} .

$N_V$	Vel.	Cell #	#nodes	$N_{\text{cycle}}$	Time(s)	$T_{\text{cycle}}$	$T_{\text{cell}}$	$T_{\text{cell/node}}$
	$[-15,15]$		15 (300 cores)		$6815.39$	$51.2$	${\mathrm{2.60\; 10}}^{-6}$	${\mathrm{3.91\; 10}}^{-5}$
			18 (360 cores)		$5479.12$	$41.2$	${\mathrm{2.09\; 10}}^{-6}$	${\mathrm{3.77\; 10}}^{-5}$
			27 (540 cores)		$3433.45$	$25.8$	${\mathrm{1.31\; 10}}^{-6}$	${\mathrm{3.54\; 10}}^{-5}$
${16}^3$		${270}^3\times {16}^3$	30 (600 cores)	$133$	$3078.33$	$23.1$	${\mathrm{1.18\; 10}}^{-6}$	${\mathrm{3.53\; 10}}^{-5}$
		$=80.6\times {10}^9$	45 (900 cores)		$2089.94$	$15.7$	${\mathrm{7.98\; 10}}^{-7}$	${\mathrm{3.59\; 10}}^{-5}$
			54 (1080 cores)		$1758$	$13.2$	${\mathrm{6.72\; 10}}^{-7}$	${\mathrm{3.63\; 10}}^{-5}$
			90 (1800 cores)		$1112.93$	$8.37$	${\mathrm{4.25\; 10}}^{-7}$	${\mathrm{3.83\; 10}}^{-5}$

Table 3. Table 3: Performance tests on 192 3 × 16 3 superscript 192 3 superscript 16 3 192^{3}\times 16^{3} mesh for BGK. Time per cycle is obtained by T cycle = T / N cycle subscript 𝑇 cycle 𝑇 subscript 𝑁 cycle T_{\text{cycle}}=T/N_{\text{cycle}} , time per cycle per cell by T cell = T cycle / N c subscript 𝑇 cell subscript 𝑇 cycle subscript 𝑁 𝑐 T_{\text{cell}}=T_{\text{cycle}}/N_{c} and time per cycle per node by T cycle/node = N s ​ T cell subscript 𝑇 cycle/node subscript 𝑁 𝑠 subscript 𝑇 cell T_{\text{cycle/node}}=N_{s}T_{\text{cell}} .

$N_V$	Vel.	Cell #	#nodes	$N_{\text{cycle}}$	Time(s)	$T_{\text{cycle}}$	$T_{\text{cell}}$	$T_{\text{cycle/node}}$
	$[-15,15]$		4 (80 cores)		$6696.82$	$70.5$	${\mathrm{9.96\; 10}}^{-6}$	${\mathrm{3.98\; 10}}^{-5}$
			6 (120 cores)		$4322.4$	$45.5$	${\mathrm{6.43\; 10}}^{-6}$	${\mathrm{3.86\; 10}}^{-5}$
			8 (160 cores)		$3012.17$	$31.7$	${\mathrm{4.48\; 10}}^{-6}$	${\mathrm{3.58\; 10}}^{-5}$
			12 (240 cores)		$1954.71$	$20.6$	${\mathrm{2.91\; 10}}^{-6}$	${\mathrm{3.49\; 10}}^{-5}$
${16}^3$		${192}^3\times {16}^3$	16 (320 cores)	$95$	$1482.56$	$15.6$	${\mathrm{2.20\; 10}}^{-6}$	${\mathrm{3.53\; 10}}^{-5}$
		$=29\times {10}^9$	24 (480 cores)		$1000.37$	$10.5$	${\mathrm{1.49\; 10}}^{-6}$	${\mathrm{3.57\; 10}}^{-5}$
			32 (640 cores)		$768.273$	$8.09$	${\mathrm{1.14\; 10}}^{-6}$	${\mathrm{3.66\; 10}}^{-5}$
			48 (960 cores)		$529.223$	$5.57$	${\mathrm{7.87\; 10}}^{-7}$	${\mathrm{3.78\; 10}}^{-5}$
			64 (1280 cores)		$409.52$	$4.31$	${\mathrm{6.09\; 10}}^{-7}$	${\mathrm{3.90\; 10}}^{-5}$

Table 4. Table 4: Performance tests on 128 3 × 16 3 superscript 128 3 superscript 16 3 128^{3}\times 16^{3} mesh for BGK. Time per cycle is obtained by T cycle = T / N cycle subscript 𝑇 cycle 𝑇 subscript 𝑁 cycle T_{\text{cycle}}=T/N_{\text{cycle}} , time per cycle per cell by T cell = T cycle / N c subscript 𝑇 cell subscript 𝑇 cycle subscript 𝑁 𝑐 T_{\text{cell}}=T_{\text{cycle}}/N_{c} and time per cycle per node by T cycle/node = N s ​ T cell subscript 𝑇 cycle/node subscript 𝑁 𝑠 subscript 𝑇 cell T_{\text{cycle/node}}=N_{s}T_{\text{cell}} .

$N_V$	Vel.	Cell #	#nodes	$N_{\text{cycle}}$	Time(s)	$T_{\text{cycle}}$	$T_{\text{cell}}$	$T_{\text{cell/node}}$
	$[-15,15]$		2 (40 cores)		$2499.49$	$39.1$	${\mathrm{1.86\; 10}}^{-5}$	${\mathrm{3.72\; 10}}^{-5}$
			4 (80 cores)		$1167.2$	$18.2$	${\mathrm{8.70\; 10}}^{-6}$	${\mathrm{3.48\; 10}}^{-5}$
${16}^3$		${128}^3\times {16}^3$	8 (160 cores)	$64$	$592.256$	$9.25$	${\mathrm{4.41\; 10}}^{-6}$	${\mathrm{3.53\; 10}}^{-5}$
		$=8.6\times {10}^9$	16 (320 cores)		$306.597$	$4.79$	${\mathrm{2.28\; 10}}^{-6}$	${\mathrm{3.65\; 10}}^{-5}$
			32 (640 cores)		$160.077$	$2.50$	${\mathrm{1.19\; 10}}^{-6}$	${\mathrm{3.82\; 10}}^{-5}$
			64 (1280 cores)		$88.1792$	$1.38$	${\mathrm{6.57\; 10}}^{-7}$	${\mathrm{4.20\; 10}}^{-5}$

Table 5. Table 5: Performance tests on 64 3 × 32 3 superscript 64 3 superscript 32 3 64^{3}\times 32^{3} mesh for BGK. Time per cycle is obtained by T cycle = T / N cycle subscript 𝑇 cycle 𝑇 subscript 𝑁 cycle T_{\text{cycle}}=T/N_{\text{cycle}} , time per cycle per cell by T cell = T cycle / N c subscript 𝑇 cell subscript 𝑇 cycle subscript 𝑁 𝑐 T_{\text{cell}}=T_{\text{cycle}}/N_{c} and time per cycle per node by T cycle/node = N s ​ T cell subscript 𝑇 cycle/node subscript 𝑁 𝑠 subscript 𝑇 cell T_{\text{cycle/node}}=N_{s}T_{\text{cell}} .

$N_V$	Vel.	Cell #	#nodes	$N_{\text{cycle}}$	Time(s)	$T_{\text{cycle}}$	$T_{\text{cell}}$	$T_{\text{cell/node}}$
	$[-15,15]$		2 (40 cores)		$1250.51$	$37.9$	${\mathrm{1.45\; 10}}^{-4}$	${\mathrm{2.89\; 10}}^{-4}$
			4 (80 cores)		$635.297$	$19.3$	${\mathrm{7.34\; 10}}^{-5}$	${\mathrm{2.94\; 10}}^{-4}$
${32}^3$		${64}^3\times {32}^3$	8 (160 cores)	$33$	$328.987$	$9.97$	${\mathrm{3.80\; 10}}^{-5}$	${\mathrm{3.04\; 10}}^{-4}$
		$=8.6\times {10}^9$	16 (320 cores)		$172.168$	$5.22$	${\mathrm{1.99\; 10}}^{-5}$	${\mathrm{3.18\; 10}}^{-4}$
			32 (640 cores)		$98.2927$	$2.98$	${\mathrm{1.14\; 10}}^{-5}$	${\mathrm{3.64\; 10}}^{-4}$
			64 (1280 cores)		$58.6773$	$1.78$	${\mathrm{6.78\; 10}}^{-6}$	${\mathrm{4.34\; 10}}^{-4}$

Table 6. Table 6: Profiling of the average cost for each routine for the BGK collision operator on 64 3 × 32 3 superscript 64 3 superscript 32 3 64^{3}\times 32^{3} mesh for 32 time steps and for 270 3 × 16 3 superscript 270 3 superscript 16 3 270^{3}\times 16^{3} mesh for 133 time steps.

		Cycle	CPU	Main routines	Cost CPU	Cost
			(s)		vs total (s)	vs total (%)
${64}^3\times {32}^3$	2 nodes	32	1374	Transport	0.0028	${\mathrm{2\; 10}}^{-4}$%
				ToConservative	269.123	19.59%
				ToPrimitive	0.01	${\mathrm{7\; 10}}^{-4}$%
				Collision	904.166	65.8%
				Communication	200.718	14.61%
				=	1374.02	100%
${64}^3\times {32}^3$	64 nodes	32	56	Transport	0.0027	${\mathrm{4.9\; 10}}^{-3}$%
				ToConservative	18.06	32.34%
				ToPrimitive	0.0008	${\mathrm{1.4\; 10}}^{-3}$%
				Collision	28.113	50.34%
				Communication	9.67	17.31%
				=	55.85	100%
${270}^3\times {16}^3$	90 nodes	133	1038	Transport	0.0041	${\mathrm{4\; 10}}^{-4}$%
				ToConservative	234.78	22.62%
				ToPrimitive	0.053	${\mathrm{5\; 10}}^{-3}$%
				Collision	751.36	72.4%
				Communication	51.59	4.97%
				=	1037.76	100%

Table 7. Table 7: Performance tests on 64 3 × 16 3 superscript 64 3 superscript 16 3 64^{3}\times 16^{3} mesh for Boltzmann collision kernel. Time per cycle is obtained by T cycle = T / N cycle subscript 𝑇 cycle 𝑇 subscript 𝑁 cycle T_{\text{cycle}}=T/N_{\text{cycle}} , time per cycle per cell by T cell = T cycle / N c subscript 𝑇 cell subscript 𝑇 cycle subscript 𝑁 𝑐 T_{\text{cell}}=T_{\text{cycle}}/N_{c} and time per cycle per node by T cycle/node = N s ​ T cell subscript 𝑇 cycle/node subscript 𝑁 𝑠 subscript 𝑇 cell T_{\text{cycle/node}}=N_{s}T_{\text{cell}} .

$N_V$	Vel.	Cell #	#nodes	$N_{\text{cycle}}$	Time(s)	$T_{\text{cycle}}$	$T_{\text{cell}}$	$T_{\text{cell/node}}$
	$[-15,15]$		2 (40 cores)		$8043.07$	$244$	${\mathrm{9.3\; 10}}^{-4}$	${\mathrm{1.86\; 10}}^{-3}$
			4 (80 cores)		$4119.41$	$125$	${\mathrm{4.76\; 10}}^{-4}$	${\mathrm{1.90\; 10}}^{-3}$
${16}^3$		${64}^3\times {16}^3$	8 (160 cores)	$33$	$2053.01$	$62.2$	${\mathrm{2.37\; 10}}^{-4}$	${\mathrm{1.90\; 10}}^{-3}$
		$=1.07\times {10}^9$	16 (320 cores)		$1031.67$	$31.3$	${\mathrm{1.19\; 10}}^{-4}$	${\mathrm{1.91\; 10}}^{-3}$
			32 (640 cores)		$529.246$	$16$	${\mathrm{6.12\; 10}}^{-5}$	${\mathrm{1.96\; 10}}^{-3}$
			64 (1280 cores)		$264.725$	$8.02$	${\mathrm{3.06\; 10}}^{-5}$	${\mathrm{1.96\; 10}}^{-3}$

Table 8. Table 8: Performance tests on 64 3 × 32 3 superscript 64 3 superscript 32 3 64^{3}\times 32^{3} mesh for Boltzmann collision kernel. Time per cycle is obtained by T cycle = T / N cycle subscript 𝑇 cycle 𝑇 subscript 𝑁 cycle T_{\text{cycle}}=T/N_{\text{cycle}} , time per cycle per cell by T cell = T cycle / N c subscript 𝑇 cell subscript 𝑇 cycle subscript 𝑁 𝑐 T_{\text{cell}}=T_{\text{cycle}}/N_{c} and time per cycle per node by T cycle/node = N s ​ T cell subscript 𝑇 cycle/node subscript 𝑁 𝑠 subscript 𝑇 cell T_{\text{cycle/node}}=N_{s}T_{\text{cell}} .

$N_V$	Vel.	Cell #	#nodes	$N_{\text{cycle}}$	Time(s)	$T_{\text{cycle}}$	$T_{\text{cell}}$	$T_{\text{cell/node}}$
	$[-15,15]$		2 (40 cores)		$128378$	$3890$	${\mathrm{1.48\; 10}}^{-2}$	${\mathrm{2.97\; 10}}^{-2}$
			4 (80 cores)		$64233.2$	$1950$	${\mathrm{7.43\; 10}}^{-3}$	${\mathrm{2.97\; 10}}^{-2}$
${32}^3$		${64}^3\times {32}^3$	8 (160 cores)	$33$	$32506.4$	$985$	${\mathrm{3.76\; 10}}^{-3}$	${\mathrm{3.01\; 10}}^{-2}$
		$=8.6\times {10}^9$	16 (320 cores)		$16133.7$	$489$	${\mathrm{1.87\; 10}}^{-3}$	${\mathrm{2.98\; 10}}^{-2}$
			32 (640 cores)		$8194.58$	$248$	${\mathrm{9.47\; 10}}^{-4}$	${\mathrm{3.03\; 10}}^{-2}$
			64 (1280 cores)		$4088.34$	$124$	${\mathrm{4.73\; 10}}^{-4}$	${\mathrm{3.02\; 10}}^{-2}$

Table 9. Table 9: Performance tests on 128 3 × 32 3 superscript 128 3 superscript 32 3 128^{3}\times 32^{3} mesh for Boltzmann collision kernel for 16 time steps. Time per cycle is obtained by T cycle = T / N cycle subscript 𝑇 cycle 𝑇 subscript 𝑁 cycle T_{\text{cycle}}=T/N_{\text{cycle}} , time per cycle per cell by T cell = T cycle / N c subscript 𝑇 cell subscript 𝑇 cycle subscript 𝑁 𝑐 T_{\text{cell}}=T_{\text{cycle}}/N_{c} and time per cycle per node by T cycle/node = N s ​ T cell subscript 𝑇 cycle/node subscript 𝑁 𝑠 subscript 𝑇 cell T_{\text{cycle/node}}=N_{s}T_{\text{cell}} .

$N_V$	Vel.	Cell #	#nodes	$N_{\text{cycle}}$	Time(s)	$T_{\text{cycle}}$	$T_{\text{cell}}$	$T_{\text{cell/node}}$
	$[-15,15]$		16 (256 cores)		$61346.4$	$3834.2$	${\mathrm{1.83\; 10}}^{-3}$	${\mathrm{2.93\; 10}}^{-2}$
			32 (512 cores)		$31006.9$	$1937.9$	${\mathrm{9.24\; 10}}^{-4}$	${\mathrm{2.96\; 10}}^{-2}$
			64 (1024 cores)		$15448.8$	$965.6$	${\mathrm{4.60\; 10}}^{-4}$	${\mathrm{2.95\; 10}}^{-2}$
${32}^3$		${128}^3\times {32}^3$	128 (2048 cores)	$16$	$7852.4$	$490.8$	${\mathrm{2.34\; 10}}^{-4}$	${\mathrm{3.00\; 10}}^{-2}$
		$=69\times {10}^9$	256 (4096 cores)		$3924.8$	$245.3$	${\mathrm{1.17\; 10}}^{-4}$	${\mathrm{2.99\; 10}}^{-2}$
			512 (8192 cores)		$1963.7$	$122.7$	${\mathrm{5.85\; 10}}^{-5}$	${\mathrm{3.00\; 10}}^{-2}$
			1024 (16384 cores)		$987.37$	$61.71$	${\mathrm{2.94\; 10}}^{-5}$	${\mathrm{3.01\; 10}}^{-2}$

Table 10. Table 10: Profiling of the average cost for each routine for the Boltzmann collision kernel on 128 3 × 32 3 superscript 128 3 superscript 32 3 128^{3}\times 32^{3} mesh for 16 time steps.

	Cycle	CPU	Main routines	Cost CPU	Cost
		(s)		vs total (s)	vs total (%)
16 nodes	16	61346	Transport	0.0032	${\mathrm{5\; 10}}^{-6}$%
			ToConservative	175.2	0.29%
			ToPrimitive	0.0056	${\mathrm{9\; 10}}^{-6}$%
			Collision	60418.6	98.5%
			Communication	752.59	1.23%
			=	61346.4	100%
1024 nodes	16	987	Transport	0.0043	${\mathrm{4.3\; 10}}^{-3}$%
			ToConservative	6.16	0.62%
			ToPrimitive	0.00028	${\mathrm{2.9\; 10}}^{-6}$%
			Collision	944.222	95.63%
			Communication	36.98	3.75%
			=	987.365	100%

Equations

$${\partial }_{t}f+v\cdot {\nabla }_{x}f=\mathcal{Q}(f,f),$$

$$\displaystyle U=\left(\begin{array}[]{c}\rho\\ \rho u\\ E\end{array}\right)=\displaystyle\int_{\mathbb{R}^{3}}\phi(v)fdv,$$

$${\int }_{{\mathbb{R}}^3}\varphi (v)\mathcal{Q}(f,f)dv=0$$

$$\frac{\partial }{\partial t}{\int }_{{\mathbb{R}}^3}f\varphi (v)dv+{\int }_{{\mathbb{R}}^3}v\cdot {\nabla }_{x}f\varphi (v)dv=0.$$

$$M[f]=\frac{\rho }{(2\pi T{)}^{3/2}}{e}^{-\frac{(u-v{)}^2}{2T}},$$

$$\frac{3}{2}\rho T=E-\frac{1}{2}\rho |u{|}^2.$$

$$\frac{\partial }{\partial t}\rho +{\nabla }_x\cdot (\rho u)=0,$$

$$\frac{\partial (\rho u)}{\partial t}+{\nabla }_x\cdot (\rho u\otimes u+pI)=0,$$

$$\frac{\partial E}{\partial t}+{\nabla }_x\cdot \left((E+p)u\right)=0,$$

$${\mathcal{Q}}_{BGK}(f,f)=\nu (M[f]-f),$$

$${\mathcal{Q}}_B(f,f)={\int }_{{\mathbb{R}}^3}{\int }_{S^2}B(|v-v_{\star }|,\theta )\left(f(v^{\prime })f(v_{\star }^{\prime })-f(v)f(v_{\star })\right)d\omega d{v}_{\star },$$

$$v^{\prime }=\frac{1}{2}(v+v_{\star }+|v-v_{\star }|\omega ),v_{\star }^{\prime }=\frac{1}{2}(v+v_{\star }-|v-v_{\star }|\omega ),$$

$$B(|v-v_{\star }|,\theta )=|v-v_{\star }|\sigma (|v-v_{\star }|,\theta ),$$

$$\sigma (|v-v_{\star }|,\theta )=b_{\alpha }(\theta )|v-v_{\star }{|}^{\alpha -1}$$

$${\tilde{f}}_k(x,t)\approx f(x,v_k,t),$$

$${\partial }_t{\tilde{f}}_k+v_k\cdot {\nabla }_x{\tilde{f}}_k={\mathcal{Q}}_k(\tilde{f},\tilde{f}),$$

$$⟶$$

$${\bar{f}}_k^{\star ,n+1}={\bar{f}}_k^n(x-v_k\Delta t).$$

$${\partial }_t{f}_{j,k}={\mathcal{Q}}_k(f_{j,\cdot },f_{j,\cdot }),$$

$$f_{j,k}^{n+1}=f_{j,k}^{\star ,n+1}+\Delta t{\mathcal{Q}}_k(f_{j,\cdot }^{\star ,n+1},f_{j,\cdot }^{\star ,n+1}),$$

$${\bar{\mathcal{Q}}}_k^{n+1}(x)={\mathcal{Q}}_k(f_{j,\cdot }^{\star ,n+1},f_{j,\cdot }^{\star ,n+1}),\forall x\text{such that}{\bar{f}}_k^{\star ,n+1}(x)=f_{j,k}^{\star ,n+1},$$

$${\bar{f}}_k^{n+1}(x)={\bar{f}}_k^{\star ,n+1}(x)+\Delta t{\bar{\mathcal{Q}}}_k^{n+1}(x).$$

$$f(x,v,t)=\sum _{j,k=1}^{N_s,N_v}{\mathfrak{m}}_{j,k}(t)\delta (x-x_{j,k}(t))\delta (v-v_{j,k}(t)),v_{j,k}(t)=v_k,$$

$$x_{j,k}(t+\Delta t)=x_{j,k}(t)+v_{j,k}(t)\Delta t.$$

$${\mathfrak{m}}_{j,k}(t+\Delta t)={\mathfrak{m}}_{j,k}(t)+\Delta t{Q}_k(v_{j,\cdot }),$$

$$U_j^n=\sum _{k=0}^{N_v}\varphi (v_{j,k}){\mathfrak{m}}_{j,k}^n(\Delta v{)}^3.$$

$$U_j^{n+1}=U_j^n+\sum _{k,x_{j,k}^{n+1}\in {\Omega }_{j+\delta },x_{j,k}^n\in {\Omega }_j}({\mathfrak{m}}_{j-\delta ,k}^n-{\mathfrak{m}}_{j,k}^n)\varphi (v_{j,k})(\Delta v{)}^3$$

$${\mathcal{Q}}_B(f,f)={\int }_{B_0(2R)}{\int }_{S^2}B(|g|,\theta )\left(f(v^{\prime })f(v_1^{\prime })-f(v)f(v-g)\right)d\omega dg$$

$$f_N(v)=\sum _{k=-N/2}^{N/2}{\hat{f}}_k{e}^{ik\cdot v},$$

$${\hat{f}}_k=\frac{1}{(2\pi {)}^3}{\int }_{[-\pi ,\pi {]}^3}f(v)e^{-ik\cdot v}dv.$$

$${\int }_{[-\pi ,\pi {]}^3}\left(\frac{\partial f_N}{\partial t}-Q_B^R(f_N)\right)e^{-ik\cdot v}dv=0.$$

Full text

Keywords: Boltzmann equation, kinetic equations, semi-Lagrangian schemes, spectral schemes, 3D/3D, MPI

1 Introduction

Kinetic equations provide a statistical description of non equilibrium particle gases. The evolution of the system is described by a ballistic motion of particles interacting only by a two-body collisions [10, 20]. The Boltzmann model derived originally in 1870s for rarefied gases that are far from thermodynamic equilibrium is nowadays used in variety of applications: ranging from plasma physics to astrophysics, quantum physics, biology and social science. In the Boltzmann description, the state of the system is described by a distribution function defined in seven independent dimensions: the physical space, the velocity space and the time. Moreover, the interaction term requires multiple integrals over velocity space to be evaluated at every space point and for every time step of the numerical method [22, 40]. This makes the kinetic theory very challenging from numerical view point.

There are two major strategies to approach numerically the Boltzmann equation. The first is to apply probabilistic methods such as Direct Simulation Monte Carlo (DSMC) [4, 7, 8, 37]. The second is to choose a deterministic scheme such as finite volume or spectral methods [20, 23, 25, 35, 38]. The probabilistic approach more efficient in terms of computational time. It is however only low order with slow convergence rate.

In this work we choose a semi-Lagrangian approach [11, 12, 23, 24, 32, 42] applied to the transport part of the Boltzmann equation coupled with spectral methods to solve the collision operator [6, 9, 21, 22, 24, 30, 31, 29, 39, 41, 46, 47, 48, 49]. In particular, we consider a Fast Kinetic Scheme (FKS) developed originally for the the Bhatnagar-Gross-Krook (BGK) operator [3] in [13, 14, 17]. The FKS applies the Discrete Velocity Model (DVM) technique, where the velocity space is truncated and discretized with a set of fixed discrete velocities. As a result, the original continuous kinetic equation is replaced by a discrete set of transport equations that can be solved exactly in the semi-Lagrangian framework at practically no cost. In the original method for the BGK operator, where the collisions are modelled as a relaxation towards the local thermodynamic equilibrium, the coupling between equations was included in the computation of the local macroscopic variables (density, momentum, temperature) used to approximate the local equilibrium state. We extend herein the FKS solver to take into account more complex collision models, such as the Boltzmann operator [4, 10]. We make use of the fast spectral method allowing to compute the collision operator in $O(N_{v}\log N_{v})$, where $N_{v}$ is a number of discrete velocity points in three dimensions [21, 22, 36].

The curse of dimensionality makes numerical simulations of the Boltzmann equation prohibitive on sequential machines even if fast numerical schemes are employed. That is why the need for efficient parallelization strategies arises. The parallel computing in the context of kinetic equations was already explored in [26, 27, 28], where the authors made use of Graphics Processing Unit (GPU) to solve the BGK equation with probabilistic methods. In [34] the authors have implemented the Boltzmann collision kernel on GPU, OpenMP and MPI algorithms were explored in [1, 2, 33, 43, 44]. The FKS parallelization on shared memory systems under OpenMP on GPU was proposed in [15] for the BGK collision kernel. In [16] a simple parallelization strategy on distributed memory systems was proposed for the FKS coupled with Boltzmann collision kernel. The goal of this article is to study in more detail the performance of the FKS on tera- and peta-scale systems. In particular, we propose an efficient hybrid MPI/OpenMP implementation of the scheme with strong scaling close to ideal on available systems.

The article is organized as follows. Section 2 introduces the kinetic equation, the Fast Kinetic Scheme and some collision operators. Section 3 discusses the particle of the FKS interpretation as particularly well suited for parallelization, Section 4 proposes a parallel algorithm and finally Section 5 gives scalability results and some profiling information.

2 Kinetic equations and Fast Kinetic Scheme

In the kinetic theory of rarefied gases, the state of the system is described by a non negative distribution function $f(x,v,t)$. This distribution function describes a density of particles moving with the velocity $v\in\mathbb{R}^{3}$ at the position $x\in\mathbb{R}^{3}$ at time $t$. The evolution of the system us governed by the six dimensional Boltzmann equation

[TABLE]

where the operator $\mathcal{Q}(f,f)$ is the collision operator and describes the effect of the collisions on the system. The macroscopic characteristics (density, momentum and energy) are obtained by integrating the distribution function multiplied by $1$, $v$ or $|v|^{2}$ over the velocity space:

[TABLE]

where the vector $\phi(v)$ is given by $(1,v,\frac{1}{2}|v|^{2})^{T}$. Typically, the collision operator conserves the macroscopic quantities of the system, i.e. the collisions preserve mass, momentum and energy. This is expressed as

[TABLE]

and hence the components of the vector $\phi(v)$ are referred to as collision invariants.

Multiplying the Boltzmann equation (1) by collision invariants and integrating over the velocity space yields a system of equations for evolution of macroscopic conservative variables

[TABLE]

This system is not closed as the second term involves higher order moments. However, when the system is at thermal equilibrium, the collision operator $\mathcal{Q}(f,f)=0$. The equilibrium is characterized by a local Maxwellian distribution

[TABLE]

where the temperature $T$ is related to the difference between the total and kinetic energy by the following relation:

[TABLE]

Replacing the distribution function $f$ in (4) by the Maxwellian distribution $M[f]$ yields a closed system — a set of Euler equations

[TABLE]

with the pressure following the ideal gas law $p=\rho T$.

The simplest collision operator providing the desired properties (conservation of collision invariant and vanishing at equilibrium) is the Bhatnagar-Gross-Krook operator [3]

[TABLE]

where the inter-particle collisions are modelled as a relaxation process towards local equilibrium. The parameter $\nu=\nu(x,t)$ defines the collision frequency.

The classical Boltzmann collision operator is a multiple integral over the whole velocity space and all possible relative angles:

[TABLE]

where $v$, $v_{\star}$ are velocities before collision, $v^{\prime}$, $v^{\prime}_{\star}$ the velocities after collision and $\theta$ the angle between $v-v_{\star}$ and $v^{\prime}-v^{\prime}_{\star}$. The post collision velocities are given by

[TABLE]

with $\omega$ being a vector on a unitary sphere $S^{2}$. The collision kernel $B$ depends only on the relative velocity before collision and the deflection angle and has the form

[TABLE]

with $\sigma$ being the scattering cross section. If the inverse $k$-th power forces between particles are considered, $\sigma$ is given by

[TABLE]

with $\alpha=(k-5)/(k-1)$. In the framework of the so called variable hard spheres model (VHS) $b_{\alpha}(\theta)$ is constant: $b_{\alpha}(\theta)=C_{\alpha}$.

2.1 Fast Kinetic Scheme

Let us now introduce the Fast Kinetic Scheme (FKS) [13, 14] for solving the Boltzmann equation (1). The FKS is a semi-Lagrangian scheme [11, 12, 25] which employs Discrete Velocity Model (DVM) [5, 35] approximation to the original problem.

Let us start by truncating a velocity space. Let us also introduce a cubic grid of $N_{v}$ equally spaced points in three dimensions. Let us assume for simplicity that the grid step $\Delta v$ is equal in every direction. Please note however that the FKS is not restricted to Cartesian grids in the velocity space. In fact the method rests unchanged even if unstructured and anisotropic velocity grids are taken into account. Let us now define an approximation of the continuous distribution function $f(x,v,t)$:

[TABLE]

that is to say, continuous $f$ is replaced by a vector $\tilde{f}$ and the following system of $N_{v}$ equations is obtained:

[TABLE]

where $\mathcal{Q}_{k}(\tilde{f},\tilde{f})$ is a suitable approximation of the collision operator for the discrete velocity point $v_{k}$. This set of equations is coupled only by the collision term.

Let us now discretize the physical space with a Cartesian grid consisting of $N_{s}$ equally spaced points with a grid step $\Delta x$ that is equal (for simplicity reasons) in all three directions. Let us also introduce a time discretization with $\Delta t$ being a time step and $t^{n}=t^{0}+n\Delta t$ for any $n\geq 0$.

In the FKS framework, the equation (7) is solved with a first order splitting technique. First the transport step exactly solves the left hand side of the problem, than the collision step introduces the interaction using the result of the transport step as a starting point:

[TABLE]

Please note that higher order splitting techniques may also be considered.

Transport step

Let $f^{n}_{j,k}$ be a point-wise approximation of the distribution function at time $t^{n}$, position $x_{j}$ and velocity $v_{k}$: $f^{n}_{j,k}=f(x_{j},v_{k},t^{n})$. The main idea behind FKS is to solve the transport step (8) exactly. Let us define a piece wise constant in space approximation $\bar{f}^{n}_{k}(x)$ of the function $\tilde{f}_{k}(x,t^{n})$ such that $\bar{f}^{n}_{k}(x)=f(x_{j},v_{k},t^{n})$ if $x\in[x_{j-1/2},x_{j+1/2}]=\Omega_{j}$ belongs to the space cell centered on $x_{j}$. The exact solution to the transport step at time $t^{n}$ is therefore given by

[TABLE]

The function $\bar{f}^{n}_{k}$ is advected with a velocity $v_{k}$ during a time step $\Delta t$. The discontinuities of $\bar{f}^{\star,n+1}_{k}$ does coincide now with space cell boundaries after the transport step.

Collision step

During the collision step the amplitude of the distribution function $\bar{f}$ is modified. The collision operator is solved locally on the space grid points and than extended to the whole domain $\Omega$. The following equations (ordinary differential or integro-differential) are solved:

[TABLE]

where $f_{j,k}=f(x_{j},v_{k},t)$ for all space and velocity grid points $j=1,\ldots,N_{s}$ and $k=1,\ldots,N_{v}$ and $f_{j,\cdot}$ is a vector representing the distribution function at the space cell $j$ composed of $f_{j,k}$. The initial data for this equation are provided by the transport step performed before. The time discretization chosen in this work is the first order explicit Euler scheme

[TABLE]

where $f_{j,k}^{\star,n+1}=\bar{f}_{k}^{\star,n+1}(x_{j})$ is the value of transported distribution function at grid point $x_{j}$ and $f_{j,\cdot}^{\star,n+1}$ is a vector composed of $f_{j,k}^{\star,n+1}$. Please note that other type of time integrators can be successfully implemented instead of this forward scheme. In particular the special care must be taken in the stiff limit, please refer to [18, 19].

Equation (9) furnishes a modified value of the distribution function at grid points $x_{j}$ for velocity points $v_{k}$ at time $t^{n+1}$. In order to obtain the value of $f$ at every point of the domain a new piecewise constant function $\bar{\mathcal{Q}}_{k}$ is defined for every discrete velocity $v_{k}$:

[TABLE]

that is to say, the collision operator at every point of $\Omega$ is approximated by a piecewise constant function with discontinuities located at the same points as the piecewise constant function that approximates the distribution function after the transport step. Thanks to this assumption, the spatially reconstructed distribution function after the collision step reads

[TABLE]

This completes the description of the Fast Kinetic Scheme.

2.2 Collision operator

Let us now focus on the details related to the collision kernel.

2.2.1 BGK approximation

If particle interaction is modeled by relaxation towards local equilibrium, the collision term $\mathcal{Q}_{k}(f_{j,\cdot}^{\star,n+1},f_{j,\cdot}^{\star,n+1})$ becomes $\nu({\cal E}_{j,k}-f_{j,k})$, where ${\cal E}_{j,k}$ is a suitable approximation of the Maxwell distribution for the velocity $v_{k}$ at the grid point $x_{j}$. As the Maxwellian distribution depends on the macroscopic characteristics of the system that are unchanged during the relaxation step (since they are collision invariants), the relaxation step (9) becomes completely decoupled. In particular, $\mathcal{Q}_{k}$ depends only on one velocity point $v_{k}$ and not on the others.

2.2.2 Boltzmann operator

If the Boltzmann operator is considered, the collision operator $\mathcal{Q}_{k}(f_{j,\cdot}^{\star,n+1},f_{j,\cdot}^{\star,n+1})$ involves integration over whole velocity space for every point $x_{j}$ of the space grid. The relaxation step is solved by means of Fast Spectral Scheme presented in the Appendix A and requires multiple Fourier transforms to be computed at every time step and at every space cell.

3 Implementation

Let us switch to a particle interpretation of the FKS. Every point of the velocity grid represents a particle moving with velocity $v_{k}$. Every space cell $\Omega_{j}$ centered on the space grid point $x_{j}$ contains exactly the same set of particles at exactly the same relative positions. Therefore one needs to store the particle position and velocity only in one generic cell and not in the whole domain. This reduces the memory requirements seven times: only mass of the particles is stored for every point of the $6D$ grid, three components of particle position and velocity vectors are only required for the generic reference cell. The distribution function is related to particle masses by

[TABLE]

where $x_{j,k}(t)$ is particle position, $v_{j,k}(t)$ is its velocity and $\mathfrak{m}_{j,k}(t)$ particle mass. In the FKS the particle velocity is unchanged and the position is altered during the transport step:

[TABLE]

The transport step moves the particles in the reference cell. The motion of particles in the remaining cells is identical. If a given particle escapes the generic cell, another one with the same velocity (but different mass) enters the cell from the opposite side.

The collision step modifies the particle masses in every space cell:

[TABLE]

where $Q_{k}(v_{j,\cdot})$ is the approximation of the collision operator in the center of the cell by the means of the fast spectral method presented above.

The macroscopic variables at time $t^{n}$ are defined on the space grid only and are computed as a sum over particles in the given cell $j$:

[TABLE]

As the collision step does not change the macroscopic conservative variables, they can be efficiently computed at time $t^{n+1}$ after the transport step by adding the contribution from the particles leaving and entering the given cell $j$ to the values at the previous time step. If a particle $(j,k)$ is transported to the cell $j+\delta$ during the transport step, there is a sister particle entering the cell $j$ from $j-\delta$. A suitable contribution has to be added to from the conservative variables in the cell $j$ :

[TABLE]

The most expensive part in the algorithm is the collision operator. Even in the case of the relatively simple BGK approximation the computation of the relaxation term takes $90\%$ of the computational time on serial machines [15]. The cost of the Boltzmann integral is substantially greater, even if the Fast Spectral Method is employed. Indeed, a number of FFTs must be performed for every space cell and for every discrete angle in order to compute convolutions. If 16 discrete angles are considered ($A_{1}=A_{2}=4$), this number equals 96 and even if the size of those transforms is relatively small, this represents the main computational burden. Evaluation of the Boltzmann operator represents more than $99\%$ of the computational time on serial machines. Fortunately the collision operator is in some sense decoupled from the FKS framework: it can be implemented independently of the FKS. This suggest a following strategy for the parallelization on distributed memory systems. On the upper level, the FKS is parallelized with MPI over available computational nodes. On the collision level, a suitable operator is implemented on the available node architecture: using the classical OpenMP type parallelism or the SIMD (Single Instruction, Multiple Data) programming model on GPUs or on the Intel Many Integrated Core (MIC) hardware. The implementation details of the collision operator does not influence the MPI scalability of the algorithm.

4 Fast Kinetic Scheme on parallel machines

In this section we propose a simple parallelization strategy for distributed memory systems. There are two possible approaches: the first one is to decompose velocity space and distribute it over computational nodes keeping all spatial degrees of freedom at every node. This approach is very similar to the strategy employed in the GPU algorithm for the BGK collision kernel (cf. [15] ), where for each velocity point from the velocity mesh a relaxation term was computed in parallel for all space cells at once on a GPU device. It was also chosen in [43] for the MPI implementation. Every computational node performs computations of a relaxation for a subset of velocity grid. Then the partial moments are evaluated. The total moments are obtained from gathering all contributions from all computational nodes. This approach is well suited for collision kernels that are local, e.g. for the BGK approximation, where the collision computed for a given position in physical and velocity space depends only on the distribution function at the same position and on total moments. As every node contains all spatial degrees of freedom, no particle escapes given computational node and no particle mass is exchanged with neighbouring nodes. The MPI communication is limited only to the partial moments. Another advantage is that even if a complicated domain is considered i.e. containing perforations, no complicated domain decomposition or load balancing techniques are required to ensure equal workloads across computational nodes. However, any collision kernel that is non-local would generate huge amount of communication between all processors. Boltzmann collision operator involve a double integral over velocity space meaning that at every iteration every node must have access to the whole velocity space. This approach is therefore not well suited for the Boltzmann operator.

The second possibility, adopted herein, is to distribute spatial degrees of freedom over computational nodes, keeping on every mode a complete velocity space. Since the update of conservative variables (density, momentum and energy) requires an exchange of particle mass with neighbouring spatial cells (and does not involve any summation over whole physical space), the internodal memory transfer is limited to merely cells located on a boundary of a subdomain. Moreover, the information is exchanged with one node only and not with every node reserved for the computation. Comparison of the two approaches can be found in [44].

The spatial domain is decomposed into equally sized non-overlapping cuboids, pencils or slabs with ghost layers. Depending on the choice, every node has to communicate with 2 (for slabs), 8 (for pencils) or 26 (for cuboids) neighboring nodes. Cuboids usually minimize the size of ghost layers but have the biggest MPI overhead as they require more calls to MPI in order to communicate with all neighbors. Thus, the cuboid domain decomposition strategy is not necessarily optimal and better results can be sometimes obtained when using pencils or slabs.

Moreover, the performance of the method can be improved by OpenMP or SIMD parallelization applied to the loops over velocity space at each node. This kind of parallelization for shared memory systems was already proposed in [15]. Let us now, for the sake of completeness, repeat the sequential and parallel algorithms for shared memory systems.

4.1 Sequential algorithms for FKS

We consider particle positions $\bm{X}_{p}^{n}$ and masses $\mathfrak{m}_{j,p}^{n}$ known at time $t^{n}$ as well as conservative $F_{j}^{n}$ and primitives variables $(\rho,\bm{U},T)_{j}^{n}$. The algorithm reads:

Transport of particles. Displace $N_{v}$ particles with (8), produce a list of $N_{\text{out}}$ particles escaping the generic cell and store the $\delta$s determining the destination and provenance of associated sister particles. 2. 2.

Update conservative variables $U_{j}^{n+1}$ using (10) and the results from the transport step. 3. 3.

Relaxation step. Compute masses of $N_{v}$ particles with a collision kernel of choice, store them in an array of the size $N_{v}\times N_{s}$.

4.2 Classical parallel architecture: Open-MP

The modified algorithm reads:

Relaxation step. Divide the number of spatial cells by the number of processors. Compute in parallel the masses of $N_{v}$ particles with a collision kernel of choice, parallelization is performed on the loop over the number of mesh points in the physical space. This computation is local on the space mesh. 2. 2.

Transport of particles. Move in parallel $N_{v}$ particles with (8). This step is done in only one space cell. The motion of particles in the other cells is the same. 3. 3.

Update conservative variables. Test in a parallel loop over the number of mesh points in the physical space if a particle has escaped from the generic cell. If so, add a contribution to $U_{j}^{n+1}$ using (10) for every space cell. Update the particle position and exchange particle mass with the associated sister particle.

4.3 Graphic Processing Unit (GPU) architecture: CUDA

This parallelization design can be summarized in the following algorithm.

Copy from CPU to GPU. Copy to the GPU memory all primitive and conservative variables. 2. 2.

Loop over $N_{v}$ particles

1. (a)

Relaxation step Compute relaxed masses of particles for every space cell using CUDA. Store the result on GPU. 2. (b)

Transport step Move every $N_{v}$ particle and test if it has escaped the generic cell. If so, store the provenance cell of the sister particle. 3. (c)

Update conservative variables. If the particle has escaped the generic cell, add contribution to conservative variables. Reassign its mass and position with the ones of the incoming sister particle. 4. (d)

Copy from GPU to CPU. Copy the resulting mass array from the GPU memory to the CPU memory. 3. 3.

Copy from GPU to CPU. Write to the CPU memory the updated conservative and primitive variables.

4.4 MPI version of FKS

The parallel algorithm is straightforward

Initialization Divide the computational domain into $N_{M\!P\!I}=N_{M\!P\!I_{x}}\times N_{M\!P\!I_{y}}\times N_{M\!P\!I_{z}}$ equally sized cuboids. Allocate memory on each computational node: arrays of the size $N_{s}/N_{M\!P\!I}$ for storing the conservative and primitive variables relative to a given subdomain and an array of the size $N_{v}\times(N_{s}/N_{M\!P\!I}+N_{\text{ghosts}})$ for storing masses relative to a given subdomain with additional ghost layers containing masses of particles in the adjacent space cells. 2. 2.

Time iterations For every computational node:

1. (a)

Relaxation step performed in parallel on GPU or with OpenMP for every particle in a given subdomain. 2. (b)

Transport of particles. Move in parallel $N_{v}$ particles with (8). This step is done in only one space cell in each subdomain. The motion of particles in the other cells is the same. 3. (c)

Communication If a particle is escaping from a given subdomain, broadcast its mass to suitable computational node. 4. (d)

Update conservative variables. If the particle has escaped the generic cell, add contribution to conservative variables. Reassign its mass and position with the ones of the incoming sister particle. For particles located on the boundary of the cuboid and escaping the subdomain use the values stored in the ghost cells in the previous step. 3. 3.

Finalization Free memory and close MPI communication.

The communication is limited only to neighboring subdomains. The amount of data to be exchanged depends only on the local mesh size and chosen MPI topology (slabs / pencils / cuboids). It does not depend on number of computational nodes employed. Weak scaling is therefore evident.

5 Performance tests

Numerical tests for the BGK collision kernel were performed on the EOS supercomputer at CALMIP, Toulouse. The supercomputer is equipped with 612 computational nodes, each of them containing two Intel${}^{\text{\textregistered}}$ Ivy Bridge $2.8$GHz 10 core CPUs and 64 GB of RAM. Each CPU was equipped with $25$MB of cache memory. The code was compiled with gcc-5.3.0 and executed on 2 to 90 computational nodes. That is to say, on 40 to 1800 computational cores in parallel. The tests for Boltzmann collision kernel were performed on the EOS supercomputer ($N=64^{3}$ meshes) and on the thin nodes of GENCI-TGCC supercomputer Curie for $N=128^{3}$. The machine is equipped with 5040 B510 Bullx nodes (called thin nodes), each containing two Intel${}^{\text{\textregistered}}$ Sandy Bridge $2.7$GHz 8 core CPUs (20MB of cache memory) and 64GB of RAM. In the case of the Boltzmann collision operator the Fast Fourier Transforms were computed by means of the fftw library, version 3.3.4. The code was executed with 20 OpenMP threads per node on EOS and 16 OpenMP threads per node on Curie.

The performance of the parallel algorithm was tested on the 3D Sod test case. The problem description is a 3D explosion problem [45], the initial state being given by the well known Sod shock tube problem. Let us consider a cubic domain of size $[0,2]^{3}$. Left and right states of the 1D Sod problem are given by a density $\rho_{L}=1$, mean velocity $\bm{U}_{L}=\bm{0}$ and temperature $\theta_{L}=5$, while $\rho_{R}=0.125$, $\bm{U}_{R}=\bm{0}$, $\theta_{R}=4$. The gas is initially in thermodynamic equilibrium. The left state is set for any cell inside a ball centered in $(1,1,1)$ and of radius $0.2$. The right state is set elsewhere. The computations are stopped at final time $t_{\text{final}}=0.07$ for the BGK operator and after $16$ iterations for the Boltzmann operator. This number of times steps ensured that every discrete particle has changed the physical cell at least once. We consider the case in which $\tau=10^{-1}$, i.e. far from the fluid limit. We are not interested in the convergence of the numerical solution but in the parallel efficiency only. Some convergence results and more interesting (from physical point of view) test cases are presented in [16].

The performance tests were run for the BGK and 3D Boltzmann collision kernels. Computation of the 3D Boltzmann collision operator is much more time consuming compared to the BGK operator. The run time is 100 times bigger for $N_{v}=32^{3}$ velocity points. The ratio increases even more for larger velocity meshes. We expect that all the communication time and MPI overhead will be hidden for the 3D Boltzmann operator with MPI efficiency close to one on the supercomputer at our disposal. When dealing with the BGK kernel it is the memory requirement that is a real bottleneck. Indeed, BGK $3D\times 3D$ simulations on medium size meshes can be efficiently performed on multi-core or GPU based shared memory systems with enough RAM at disposition with runtime that of the order of hours or days rather than weeks or months [15]. However, if fine scale simulations are required the MPI parallelization is indispensable.

5.1 BGK

Let us first discuss numerical results for the BGK collision kernel. Simulations were run on four different meshes:

• •

$N=64^{3}$ and $N_{v}=32^{3}$,

• •

$N=128^{3}$ and $N_{v}=16^{3}$,

• •

$N=192^{3}$ and $N_{v}=16^{3}$.

• •

$N=270^{3}$ and $N_{v}=16^{3}$.

The memory required for storage of the mass array was equal to $64$Gb in the first two cases, to $216$Gb in the third case and to $600$Gb in the last case. The minimal number of computational nodes required to run the Sod test case was therefore equal to $2$ for the first two cases, to $4$ in the third case and to $15$ in the last case.

Let us first compare different domain decomposition strategies for the first mesh ($N=62^{3}$ and $N_{v}=32^{3}$) to find the optimal one. The computations were run on $N_{M\!P\!I}=2$, $4$, $8$, $16$, $32$ and $64$ nodes. The domain was decomposed into slabs, pencils or cuboids. The elapsed run time is shown in Table 1. For $N_{M\!P\!I}=4$ nodes the best results are obtained for the computational domain divided into 4 slabs. This coincides with the smallest number of the ghost cells. For 8 nodes the shortest run time is still obtained for the slab configuration, even though the number of ghost cells is $28\%$ bigger than for the pencil configuration. The run time for the latter is only slightly bigger. The cuboid type domain decomposition gives the biggest run time. For 16 nodes the run time for slab and pencil like decomposition is very close to each other with relative difference of the order of $1\%$. We note that the best performance was obtained for a pencil like decomposition with $N_{M\!P\!I_{x}}=8$ and $N_{M\!P\!I_{y}}=2$. For 32 and 64 nodes the best performance was obtained again for pencil like decomposition. The cuboid type decomposition resulted in a run time comparable with the slab configuration and from $7\%$ to $12\%$ bigger than the run time for the best configuration despite the fact that the number of ghost cells was the smallest (almost $5$ time smaller than the number of ghost cells for slab configuration on 64 nodes). The weak performance of the cuboid configuration is due to a greater MPI overhead caused by communication with significantly bigger number of nodes.

The scalability test were performed with the optimal domain decomposition strategy. The run time as a function of computational cores is presented on the Figure 1 and the speedup with parallel efficiency (relative to the smallest number of nodes employed in the test) on the Figure 2. For the second mesh the scaling is almost perfect with some super-linear behaviour when passing from 2 to 4 computational nodes and some loss of performance when going from 32 to 64 for nodes. For the smallest mesh there is no super-linearity observed, the performance loss is observed when passing from 16 to 32 and then to 64 computational nodes. For the third and fourth mesh the super-linearity appears between 4 and 16 nodes employed, parallel efficiency is above $0.95$ for 2–32 nodes. The scaling for the last mesh seems to be super linear in the whole tested range with parallel efficiency above one. This super linear scaling is probably due to the CPU cache performance. For all four mesh sizes there is a linear decrease in parallel efficiency for growing number of nodes.

Tables 2–5 show the total run time $T$, time spent on a single cycle $T_{\text{cycle}}=T/N_{\text{cycle}}$, average time time spent on a single cell per cycle $T_{\text{cell}}=T_{\text{cycle}}/N_{c}$ and finally the average time spent by a single node on a single cell per cycle $T_{\text{cell/node}}=N_{M\!P\!I}T_{\text{cycle}}$ for different number of nodes employed. The performance loss is observed when the subdomain size become relatively small. This is manifested in the increase of the average time per cycle spent by one node on a single cell. This is due to the fact that the internodal communication time becomes comparable with the time spent on computations.

The size of spatial mesh allocated in every node seems also to influence the cache performance at the nodes. On Tables 2, 3 and 4 the non linear scaling is observed for small number of nodes — the average time per cell per cycle per node is decreasing when the number of nodes increases. The best performance is observed on 30 nodes for the $270^{3}$ mesh ($656\ 10^{3}$ cells per node), on 12 nodes for $192^{3}$ mesh ($590\ 10^{3}$ cells per node) and on 4 nodes on $128^{3}$ mesh ($524\ 10^{3}$ cells per node).

There is no such effect observed on Table 5. Let us now compare Tables 4 and 5. The dimension of the problem is the same. In the first case the domain was discretized with $128^{3}$ points in the physical space and with the $16^{3}$ velocity points. In the second case a spatial mesh of the size $64^{3}$ was used and a velocity mesh of the size $32^{3}$, giving a total of $8.6\ 10^{9}$ degrees of freedom. The size of the spatial mesh per node is eight times smaller in second case while the number of degrees of freedom is the same in both cases. The $N=128^{3}$ mesh seems to perform better than the $64^{3}$ one with time per cycle being smaller and with better MPI efficiency. This is not surprising as the amount of data exchanged with neighbouring nodes is at least 2 times smaller.

Table 6 shows profiling data for the $64^{3}\times 32^{3}$ mesh on 2 and 64 nodes as well as for the $270^{3}\times 16^{3}$ mesh on 90 computational nodes. For the coarse mesh the time spent on the BGK collision operator decreases from $65\%$ (on 2 nodes) to $50\%$ (on 64 nodes) while on the fine mesh the relaxation term takes more than $70\%$ of the computational time (on 90 nodes). This explains the efficiency loss for larger number of nodes: the time spent on the BGK operator is too small to hide internodal communication. Please note that better parallelization strategy can be applied for this collision model: decomposition of the velocity space (instead of physical) would require less data exchange between MPI processes in this particular case and yield therefore a better parallel efficiency. This alternative approach is however not suited for the Boltzmann collision operator.

5.2 3D Boltzmann

Let us now turn our attention to the 3D Boltzmann collision kernel. This task is much more demanding in terms of computational time as the relaxation routine involves multiple and expensive calls to the Fast Fourier Transform. The tests were run on following meshes:

• •

$N=64^{3}$ and $N_{v}=16^{3}$,

• •

$N=64^{3}$ and $N_{v}=32^{3}$.

• •

$N=128^{3}$ and $N_{v}=32^{3}$.

The tests for the $N=64^{3}$ meshes were performed on slab type domain decomposition. The results are presented on Figures 3 and 4 and on Tables 7 and 8. For the last mesh computations were performed for $16$ time steps and until the final time $t_{f}=0.07$ was reached. This clearly does not allow to present any physically interesting results but is enough to study the parallel efficiency of the herein proposed method as during this time every particle will change a space cell at least once. The tests were run on $N_{M\!P\!I}$ ranging from $16$ to $1024$. That is to say on $128$ to $2048$ processors and on $1024$ to $16384$ computational cores. The results are presented on Table 9 and on Figures 3 and 4.

A comparison of the run times obtained for the second mesh (see Figure 3) with the run time obtained for the BGK collision kernel shows that the latter is approximately 100 times faster. As the time spent on the computation of the relaxation kernel is now much more important, the communication time becomes negligible even if the size of the subdomains is small. The method shows the strong scaling that is close to ideal for all tested meshes with no efficiency loss observed (see Figure 4) and the parallel efficiency close to one. Also, the average time spent by one computational node on one cell per cycle does not depend on number of nodes employed in the computations. The method clearly enjoys the strong scaling close to ideal in the tested range.

Comparison of the run time for both meshes shows that multiplying the number of Fourier modes by 2 in each direction results in a run time multiplied by 16 instead of 8. This loss of performance is related to the computational complexity of the fast spectral solver for the Boltzmann collision kernel which is of $O(N_{v}\log N_{v})$ as well as to the cache performance. Each CPU performs 10 (8 on Curie) FFTs in parallel at the same time. The size of one transform required for the convolution computation is $32^{3}$ for the $16^{3}$ Fourier modes. This means that the memory required to store 10 vectors containing $32^{3}$ complex values in double precision format for 10 in-place transforms is $5$MB, which is less than the cache capacity. For the finer mesh the memory required is $40$MB, that is to say four times more than the cache of the CPU.

Table 10 shows the profiling data for 16 and 1024 computational nodes. Time spent on MPI communication increases from $1.23\%$ to $3.75\%$. Collision step takes above $95\%$ of the total run time and hides the parallelization costs.

6 Conclusions

In this paper we have presented an extension of the Fast Kinetic Scheme introduced in [13] to the Boltzmann collision operator by means of Fast Spectral Method and possible MPI parallelization. The obtained strong scaling is closed to ideal in the tested range (up to 1024 computational nodes) for 3D$\times$3D problems for the Boltzmann operator. In the proposed method the physical space was distributed over computational nodes by means of the MPI. Each node was supplied with a complete velocity grid. This approach has proven to be well suited for velocity non local collision operators (like Boltzmann) and performs less good for simplified and velocity local models (BGK). The reason is that the BGK collision kernel is much less resource demanding — computational time spent on particle interaction in one space cell is two orders of magnitude smaller compared to the full Boltzmann operator. Better performance for the BGK operator could be obtained by employing the following alternative: distribute velocity space over computational nodes and supply every node with a complete physical space. This approach would minimize the communication for the BGK kernel. However, the optimal algorithm for this operator is not the scope of this paper and numerical results were presented only to illustrate the computational complexity of the Boltzmann operator. The algorithm was tested on classical architectures with collision kernels parallelized locally by the means of OpenMP. We expect to maintain the similar scaling on GPU or Intel MIC based node architectures, which is the scope of future works.

Acknowledgements

This work has been supported by the french ’Agence Nationale pour la Recherche (ANR)’ in the frame of the contract “Moonrise” (ANR-11-MONU-009-01). Numerical simulations were performed using HPC resources from GENCI-TGCC (Grant 2016-AP010610045) and from CALMIP (project P1542). The author would like to acknowledge fruitful discussions with G.Dimarco, F.Filbet and R.Loubère.

Appendix A Boltzmann collision operator

The Boltzmann collision operator is solved by the means of Fast Spectral Scheme presented in this section for a selected time step $t^{n}$ and a selected grid point $x_{j}$. The same computation is repeated for other grid points and time steps. In this section we denote by $f$ the distribution function (of $v$ only) at time $t^{n}$ and point $x_{j}$: $f=f(v)=f(x_{j},v,t^{n})$.

In order to compute the Boltzmann integral (5), let us suppose that the distribution function $f$ has a compact support on the ball $B_{0}(R)$ of radius $R$ centered in the origin. It can be shown [39, 40] that the support of the collision operator $Q(f,f)$ is included in the ball $B_{0}(\sqrt{2}R)$ and

[TABLE]

with $v^{\prime},v_{1}^{\prime},v-g\in B_{0}((2+\sqrt{2})R)$. We can therefore restrict $f$ to the cube $[-T,T]^{3}$ with $T\geq(2+\sqrt{2})R$ assuming $f(v)=0$ on $[-T,T]^{3}\setminus B_{0}(R)$ and then extend it to a periodic function on $[-T,T]^{3}$. As a consequence of the periodicity of $f$ it is sufficient to take $T\geq(3+\sqrt{2})R/2$ to prevent overlapping of the regions where $f$ is different from zero [40]. In order to simplify notation let us take $T=\pi$ and $R=\lambda\pi$ with $\lambda=2/(3+\sqrt{2})$. Let $Q^{R}_{B}(f)$ denote the Boltzmann operator with cut-off. Let us perform a discrete Fourier transform on $f$ obtaining

[TABLE]

with $k$ being a multi-index $k=(k_{x},k_{y},k_{z})$ and $N=(\sqrt[3]{N_{v}},\sqrt[3]{N_{v}},\sqrt[3]{N_{v}})$ is vector containing the number of velocity discretization points (Fourier modes) in each direction. Let us impose that the residue of the collision step is orthogonal to any trigonometric polynomial of degree less or equal than $N$ in order to obtain a set of ODEs for coefficients $\hat{f}_{k}$:

[TABLE]

After some computation we obtain

[TABLE]

where $\hat{\beta}(l,m)=\hat{B}(l,m)-\hat{B}(m,m)$ are given by

[TABLE]

with

[TABLE]

Finally, the set of ODEs is obtained:

[TABLE]

supplied with the initial condition

[TABLE]

Straightforward evaluation of (11) is expensive, especially in three dimensions, the cost being of the order of $O(N^{2})$. In order to reduce the computational cost the so called Carleman representation of (5) is used:

[TABLE]

with

[TABLE]

Under the $k$-th power inter particle force assumption (6) this becomes

[TABLE]

The new quadrature formula is obtained:

[TABLE]

where ${\hat{\beta}}_{F}(l,m)=\hat{B}_{F}(l,m)-\hat{B}_{F}(m,m)$ are now given by

[TABLE]

The next step is to identify a convolution structure. The goal is to approximate ${\hat{\beta}}_{F}(l,m)$ by a sum

[TABLE]

where $A$ represents a finite number of possible collision directions. This is a discrete sum of $A$ convolutions and as a consequence the computational cost of the algorithm is of the order of $O(AN\log N)$.

This convolution can be obtained under assumption that $\tilde{B}(|x|,|y|)$ is separable:

[TABLE]

This is the case is $\alpha$ is set to one in (12): $\tilde{B}(|x|,|y|)=4C_{\alpha}$. In particular, when $\tilde{B}(|x|,|y|)=1$, this corresponds to the hard spheres model. In this framework the following quadrature formula for $\hat{B}_{F}(l,m)$ is obtained:

[TABLE]

where

[TABLE]

and the discrete angles $\theta_{p}$ and $\varphi_{q}$ are defined by

[TABLE]