A comparative study of implicit Jacobian-free Rosenbrock-Wanner, ESDIRK and BDF methods for unsteady flow simulation with high-order flux reconstruction formulations

TL;DR

This paper compares Jacobian-free Rosenbrock-Wanner, ESDIRK, and BDF2 methods for unsteady flow simulations with high-order flux reconstruction, highlighting efficiency and stability differences.

Contribution

It provides a comprehensive comparison of these methods, including implementation details and performance analysis for high-order unsteady flow simulations.

Findings

ROW and ESDIRK are more efficient than BDF2.

ROW methods can outperform ESDIRK in efficiency.

ESDIRK allows larger time steps than ROW.

Abstract

We conduct a comparative study of the Jacobian-free linearly implicit Rosenbrock-Wanner (ROW) methods, the explicit first stage, singly diagonally implicit Runge-Kutta (ESDIRK) methods, and the second-order backward differentiation formula (BDF2) for unsteady flow simulation using spatially high-order flux reconstruction/correction procedure via reconstruction (FR/CPR) formulations. The pseudo-transient continuation is employed to solve the nonlinear systems resulting from the temporal discretizations with ESDIRK and BDF2. A Jacobian-free implementation of the restarted generalized minimal residual method (GMRES) solver is employed with a low storage element-Jacobi preconditioner to solve linear systems, including those in linearly implicit ROW methods and those from linearization of the nonlinear systems in ESDIRK and BDF2 methods. We observe that all ROW and ESDIRK schemes (from second order to fourth order) are more computationally efficient than BDF2, and ROW methods can potentially be more efficient than ESDIRK methods. However, the convergence tolerance of the GMRES solver for ROW methods needs to be sufficiently tight to preserve the nominal order of accuracy. In general, ESDIRK methods allow a larger physical time step size for unsteady flow simulation than ROW methods do.

Tables

Table 1. Table 1: The convergence study for ESDIRK and BDF2 with different t ​ o ​ l r ​ e ​ l , n ​ o ​ n ​ l ​ i ​ n ​ e ​ a ​ r 𝑡 𝑜 subscript 𝑙 𝑟 𝑒 𝑙 𝑛 𝑜 𝑛 𝑙 𝑖 𝑛 𝑒 𝑎 𝑟 tol_{rel,nonlinear} when t ​ o ​ l r ​ e ​ l , l ​ i ​ n ​ e ​ a ​ r 𝑡 𝑜 subscript 𝑙 𝑟 𝑒 𝑙 𝑙 𝑖 𝑛 𝑒 𝑎 𝑟 tol_{rel,linear} is set to 10 − 1 superscript 10 1 10^{-1} .

$to{l}_{rel,nonlinear}={10}^{-2}$, $to{l}_{rel,linear}={10}^{-1}$
	$\Delta t$	$E_{L_2}\left(\rho \right)$	order	$E_{L_2}\left(u\right)$	order	CPU time(s)
ESDIRK2	T/100	$4.4097\times {10}^{-4}$		$5.7334\times {10}^{-3}$		317
	T/200	$1.1313\times {10}^{-4}$	1.97	$1.4622\times {10}^{-3}$	1.97	603
	T/400	$2.8315\times {10}^{-5}$	2.00	$3.6523\times {10}^{-4}$	2.00	1159
	T/800	$7.0714\times {10}^{-6}$	2.00	$9.1104\times {10}^{-5}$	2.00	2247
ESDIRK3	T/100	$1.5341\times {10}^{-4}$		$1.7108\times {10}^{-3}$		398
	T/200	$2.3518\times {10}^{-5}$	2.71	$2.6017\times {10}^{-4}$	2.90	734
	T/400	$3.1802\times {10}^{-6}$	2.89	$3.4016\times {10}^{-5}$	2.94	1392
	T/800	$1.0661\times {10}^{-6}$	1.58	$9.7414\times {10}^{-6}$	1.80	2478
ESDIRK4	T/100	$4.5157\times {10}^{-6}$		$3.2613\times {10}^{-5}$		425
	T/200	$7.7894\times {10}^{-7}$	2.54	$4.7845\times {10}^{-5}$	2.77	773
	T/400	$6.0568\times {10}^{-8}$	3.68	$3.2420\times {10}^{-7}$	3.88	1511
	T/800	$2.2122\times {10}^{-8}$	1.45	$1.5485\times {10}^{-7}$	1.07	2802
BDF2	T/100	$1.9170\times {10}^{-3}$		$3.0021\times {10}^{-2}$		292
	T/200	$8.1561\times {10}^{-3}$	1.23	$1.1195\times {10}^{-2}$	1.42	548
	T/400	$2.1372\times {10}^{-4}$	1.93	$2.8983\times {10}^{-3}$	1.95	1074
	T/800	$5.6824\times {10}^{-5}$	1.91	$7.3097\times {10}^{-4}$	1.99	2133
$to{l}_{rel,nonlinear}={10}^{-4}$, $to{l}_{rel,linear}={10}^{-1}$
	$\Delta t$	$E_{L_2}\left(\rho \right)$	order	$E_{L_2}\left(u\right)$	order	CPU time(s)
ESDIRK2	T/100	$4.4088\times {10}^{-4}$		$5.7335\times {10}^{-3}$		375
	T/200	$1.1278\times {10}^{-4}$	1.97	$1.4549\times {10}^{-3}$	1.98	723
	T/400	$2.8267\times {10}^{-5}$	2.00	$3.6420\times {10}^{-4}$	2.00	1323
	T/800	$7.0691\times {10}^{-6}$	2.00	$9.1040\times {10}^{-5}$	2.00	2504
ESDIRK3	T/100	$1.5167\times {10}^{-4}$		$1.6911\times {10}^{-3}$		537
	T/200	$2.3024\times {10}^{-5}$	2.72	$2.5401\times {10}^{-4}$	2.74	956
	T/400	$3.0921\times {10}^{-6}$	2.90	$3.3171\times {10}^{-5}$	2.94	1685
	T/800	$3.9376\times {10}^{-7}$	2.97	$4.1908\times {10}^{-6}$	2.98	2890
ESDIRK4	T/100	$3.4484\times {10}^{-6}$		$3.1216\times {10}^{-5}$		597
	T/200	$2.1713\times {10}^{-7}$	3.99	$1.9593\times {10}^{-6}$	4.00	1046
	T/400	$1.3573\times {10}^{-8}$	4.00	$1.2257\times {10}^{-7}$	4.00	1920
	T/800	$1.1691\times {10}^{-9}$	3.54	$7.6729\times {10}^{-9}$	3.81	3193
BDF2	T/100	$1.9126\times {10}^{-3}$		$3.0095\times {10}^{-2}$		329
	T/200	$8.4137\times {10}^{-4}$	1.18	$1.1454\times {10}^{-2}$	1.39	627
	T/400	$2.3001\times {10}^{-4}$	1.87	$2.9892\times {10}^{-3}$	1.94	1179
	T/800	$5.8417\times {10}^{-5}$	1.98	$7.5360\times {10}^{-4}$	1.99	2294
$to{l}_{rel,nonlinear}={10}^{-6}$, $to{l}_{rel,linear}={10}^{-1}$
	$\Delta t$	$E_{L_2}\left(\rho \right)$	order	$E_{L_2}\left(u\right)$	order	CPU time(s)
ESDIRK2	T/100	$4.4088\times {10}^{-4}$		$5.7335\times {10}^{-3}$		464
	T/200	$1.1278\times {10}^{-4}$	1.97	$1.4549\times {10}^{-3}$	1.98	814
	T/400	$2.8267\times {10}^{-5}$	2.00	$3.6420\times {10}^{-4}$	2.00	1503
	T/800	$7.0691\times {10}^{-6}$	2.00	$9.1041\times {10}^{-5}$	2.00	2672
ESDIRK3	T/100	$1.5165\times {10}^{-4}$		$1.6911\times {10}^{-3}$		692
	T/200	$2.3021\times {10}^{-5}$	2.72	$2.5399\times {10}^{-4}$	2.74	1138
	T/400	$3.0914\times {10}^{-6}$	2.90	$3.3165\times {10}^{-5}$	2.94	1877
	T/800	$3.9326\times {10}^{-7}$	2.97	$4.1867\times {10}^{-5}$	2.99	3369
ESDIRK4	T/100	$3.4502\times {10}^{-6}$		$3.1249\times {10}^{-5}$		780
	T/200	$2.1706\times {10}^{-7}$	3.99	$1.9593\times {10}^{-6}$	4.00	1310
	T/400	$1.3584\times {10}^{-8}$	4.00	$1.2270\times {10}^{-7}$	4.00	2156
	T/800	$8.4961\times {10}^{-10}$	4.00	$7.6728\times {10}^{-9}$	4.00	3751
BDF2	T/100	$1.9126\times {10}^{-3}$		$3.0096\times {10}^{-2}$		389
	T/200	$8.4141\times {10}^{-4}$	1.18	$1.1454\times {10}^{-2}$	1.39	708
	T/400	$2.3010\times {10}^{-4}$	1.87	$2.9890\times {10}^{-3}$	1.94	1230
	T/800	$5.8409\times {10}^{-5}$	1.98	$7.5350\times {10}^{-3}$	1.99	2463
$to{l}_{rel,nonlinear}={10}^{-8}$, $to{l}_{rel,linear}={10}^{-1}$
	$\Delta t$	$E_{L_2}\left(\rho \right)$	order	$E_{L_2}\left(u\right)$	order	CPU time(s)
ESDIRK2	T/100	$4.4088\times {10}^{-4}$		$5.7335\times {10}^{-3}$		546
	T/200	$1.1278\times {10}^{-4}$	1.97	$1.4549\times {10}^{-3}$	1.98	930
	T/400	$2.8267\times {10}^{-5}$	2.00	$3.6420\times {10}^{-4}$	2.00	1601
	T/800	$7.0691\times {10}^{-6}$	2.00	$9.1041\times {10}^{-5}$	2.00	3003
ESDIRK3	T/100	$1.5165\times {10}^{-4}$		$1.6911\times {10}^{-3}$		891
	T/200	$2.3021\times {10}^{-5}$	2.72	$2.5399\times {10}^{-4}$	2.74	1373
	T/400	$3.0913\times {10}^{-6}$	2.90	$3.3164\times {10}^{-5}$	2.94	2138
	T/800	$3.9326\times {10}^{-7}$	2.97	$4.1867\times {10}^{-5}$	2.99	3867
ESDIRK4	T/100	$3.4502\times {10}^{-6}$		$3.1249\times {10}^{-5}$		971
	T/200	$2.1706\times {10}^{-7}$	3.99	$1.9593\times {10}^{-6}$	4.00	1550
	T/400	$1.3584\times {10}^{-8}$	4.00	$1.2270\times {10}^{-7}$	4.00	2543
	T/800	$8.4976\times {10}^{-10}$	4.00	$7.6747\times {10}^{-9}$	4.00	4370
BDF2	T/100	$1.9126\times {10}^{-3}$		$3.0096\times {10}^{-2}$		455
	T/200	$8.4141\times {10}^{-4}$	1.18	$1.1454\times {10}^{-2}$	1.39	786
	T/400	$2.3010\times {10}^{-4}$	1.87	$2.9890\times {10}^{-3}$	1.94	1393
	T/800	$5.8409\times {10}^{-5}$	1.98	$7.5350\times {10}^{-3}$	1.99	2651
End of Table LABEL:fullyimplicit_nonlinear.

Table 2. Table 2: The effect of t ​ o ​ l r ​ e ​ l , l ​ i ​ n ​ e ​ a ​ r 𝑡 𝑜 subscript 𝑙 𝑟 𝑒 𝑙 𝑙 𝑖 𝑛 𝑒 𝑎 𝑟 tol_{rel,linear} on computational cost of ESDIRK.

		CPU time (second)
	$to{l}_{rel,linear}$	${10}^{-1}$	${10}^{-2}$	${10}^{-4}$	${10}^{-6}$
ESDIRK2	T/100	464	503	693	1041
	T/200	814	856	1109	1551
	T/400	1503	1563	1909	2501
	T/800	2672	2875	3334	4203
ESDIRK3	T/100	692	686	1036	1041
	T/200	1138	1128	1584	2419
	T/400	1877	1865	2580	3695
	T/800	3369	3526	4430	5908
ESDIRK4	T/100	780	875	1326	2095
	T/200	1310	1410	2133	3004
	T/400	2156	2439	3277	4844
	T/800	3751	4136	5527	7636

Table 3. Table 3: The convergence study for ROW methods with different t ​ o ​ l r ​ e ​ l , l ​ i ​ n ​ e ​ a ​ r 𝑡 𝑜 subscript 𝑙 𝑟 𝑒 𝑙 𝑙 𝑖 𝑛 𝑒 𝑎 𝑟 tol_{rel,linear} .

$to{l}_{rel,linear}={10}^{-2}$
	$\Delta t$	$E_{L_2}\left(\rho \right)$	order	$E_{L_2}\left(u\right)$	order	CPU time(s)
ROW2	T/100	Diverged	–	–	–
	T/200	$2.3615\times {10}^{-4}$		$3.3198\times {10}^{-3}$		537
	T/400	$9.9393\times {10}^{-5}$	1.25	$1.3753\times {10}^{-3}$	1.27	1052
	T/800	$5.5248\times {10}^{-5}$	0.85	$7.4627\times {10}^{-4}$	0.88	2073
ROW3	T/100	$2.6132\times {10}^{-4}$		$2.6888\times {10}^{-3}$		300
	T/200	$8.8493\times {10}^{-5}$	1.56	$8.3870\times {10}^{-4}$	1.68	564
	T/400	$3.3446\times {10}^{-5}$	1.40	$4.0554\times {10}^{-4}$	1.05	1124
	T/800	$9.1969\times {10}^{-6}$	1.86	$1.1447\times {10}^{-4}$	1.82	2191
ROW4	T/100	$3.6161\times {10}^{-5}$		$4.1926\times {10}^{-4}$		344
	T/200	$9.0325\times {10}^{-6}$	2.00	$8.3350\times {10}^{-5}$	2.33	638
	T/400	$2.7913\times {10}^{-6}$	1.69	$2.7282\times {10}^{-5}$	1.61	1194
	T/800	$3.9397\times {10}^{-6}$	-0.50	$5.1692\times {10}^{-5}$	-0.92	2282
$to{l}_{rel,linear}={10}^{-4}$
	$\Delta t$	$E_{L_2}\left(\rho \right)$	order	$E_{L_2}\left(u\right)$	order	CPU time(s)
ROW2	T/100	Diverged	–	–	–
	T/200	$8.6277\times {10}^{-5}$		$1.1163\times {10}^{-3}$		641
	T/400	$2.1800\times {10}^{-5}$	1.98	$2.8024\times {10}^{-3}$	1.99	1189
	T/800	$5.4708\times {10}^{-6}$	1.99	$7.1383\times {10}^{-5}$	1.97	2255
ROW3	T/100	$1.5467\times {10}^{-4}$		$1.7185\times {10}^{-3}$		449
	T/200	$2.3151\times {10}^{-5}$	2.74	$2.5539\times {10}^{-4}$	2.75	779
	T/400	$3.1137\times {10}^{-6}$	2.89	$3.3395\times {10}^{-5}$	2.94	1382
	T/800	$4.0729\times {10}^{-7}$	2.93	$4.2406\times {10}^{-6}$	2.98	2565
ROW4	T/100	$4.0016\times {10}^{-6}$		$3.8660\times {10}^{-5}$		487
	T/200	$3.5565\times {10}^{-7}$	3.49	$3.4725\times {10}^{-6}$	3.47	846
	T/400	$2.7298\times {10}^{-8}$	3.70	$1.6676\times {10}^{-7}$	4.38	1563
	T/800	$1.8642\times {10}^{-8}$	0.55	$2.2003\times {10}^{-7}$	-0.40	2769
$to{l}_{rel,linear}={10}^{-6}$
	$\Delta t$	$E_{L_2}\left(\rho \right)$	order	$E_{L_2}\left(u\right)$	order	CPU time(s)
ROW2	T/100	Diverged	–	–	–
	T/200	$8.6820\times {10}^{-5}$		$1.1229\times {10}^{-3}$		767
	T/400	$2.1813\times {10}^{-5}$	1.99	$2.8148\times {10}^{-4}$	2.00	1377
	T/800	$5.4628\times {10}^{-6}$	2.00	$7.0399\times {10}^{-5}$	2.00	2500
ROW3	T/100	$1.5474\times {10}^{-4}$		$1.7185\times {10}^{-3}$		690
	T/200	$2.3165\times {10}^{-5}$	2.74	$2.5534\times {10}^{-4}$	2.75	1036
	T/400	$3.1047\times {10}^{-6}$	2.90	$3.3251\times {10}^{-5}$	2.94	1712
	T/800	$4.0081\times {10}^{-7}$	2.95	$4.1768\times {10}^{-6}$	2.99	3004
ROW4	T/100	$3.9102\times {10}^{-6}$		$3.6770\times {10}^{-5}$		682
	T/200	$2.5386\times {10}^{-7}$	3.95	$2.3433\times {10}^{-6}$	3.97	1091
	T/400	$2.5819\times {10}^{-8}$	3.30	$1.3773\times {10}^{-7}$	4.09	1849
	T/800	$7.5834\times {10}^{-9}$	1.78	$1.7492\times {10}^{-8}$	2.98	3184
$to{l}_{rel,linear}={10}^{-8}$
	$\Delta t$	$E_{L_2}\left(\rho \right)$	order	$E_{L_2}\left(u\right)$	order	CPU time(s)
ROW2	T/100	Diverged	–	–	–
	T/200	$8.6824\times {10}^{-5}$		$1.1229\times {10}^{-3}$		906
	T/400	$2.1810\times {10}^{-5}$	1.99	$2.8147\times {10}^{-4}$	2.00	1557
	T/800	$5.4619\times {10}^{-6}$	2.00	$7.0398\times {10}^{-5}$	2.00	2722
ROW3	T/100	$1.5476\times {10}^{-4}$		$1.7186\times {10}^{-3}$		862
	T/200	$2.3158\times {10}^{-5}$	2.74	$2.5535\times {10}^{-4}$	2.75	1589
	T/400	$3.1015\times {10}^{-6}$	2.90	$3.3242\times {10}^{-5}$	2.94	2149
	T/800	$4.0093\times {10}^{-7}$	2.95	$4.2024\times {10}^{-6}$	2.98	3481
ROW4	T/100	$3.9055\times {10}^{-6}$		$3.6728\times {10}^{-5}$		1186
	T/200	$2.5299\times {10}^{-7}$	3.96	$2.3382\times {10}^{-6}$	3.97	1381
	T/400	$2.6107\times {10}^{-8}$	3.28	$1.4352\times {10}^{-7}$	4.03	2237
	T/800	$7.5988\times {10}^{-9}$	1.78	$1.7442\times {10}^{-8}$	3.04	3720
End of Table LABEL:linearly_implicit.

Equations

$$\frac{\partial \mathbf{q}}{\partial t}+\nabla \cdot \mathbb{F}=0,$$

$${\mathbb{F}}_{inv,j}(\mathbf{q})=\left(\begin{array}{c}\rho u_j\\ \rho u_1{u}_j+{\delta }_{1j}p\\ \rho u_2{u}_j+{\delta }_{2j}p\\ \rho u_3{u}_j+{\delta }_{3j}p\\ u_j(E+p)\end{array}\right),\text{and}{\mathbb{F}}_{vis,j}(\mathbf{q},\nabla \mathbf{q})=\left(\begin{array}{c}0\\ {\tau }_{1j}\\ {\tau }_{2j}\\ {\tau }_{3j}\\ {\sum }_{k=1}^3{u}_k{\tau }_{kj}-K_j\end{array}\right),$$

$${\tau }_{ij}=\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\frac{2}{3}\mu \sum _{k=1}^3\frac{\partial u_k}{\partial x_k}{\delta }_{ij}\text{with}i=1,2,3,$$

$$\frac{\partial \mathbf{Q}}{\partial t}+\frac{\partial \mathbf{F}}{\partial \xi }+\frac{\partial \mathbf{G}}{\partial \eta }+\frac{\partial \mathbf{H}}{\partial \zeta }=0,$$

$$\{\begin{array}{l}\mathbf{Q}=|\mathbf{J}|\mathbf{q},\\ \mathbf{F}=|\mathbf{J}|(\mathbf{f}{\xi }_x+\mathbf{g}{\xi }_y+\mathbf{h}{\xi }_z),\\ \mathbf{G}=|\mathbf{J}|(\mathbf{f}{\eta }_x+\mathbf{g}{\eta }_y+\mathbf{h}{\eta }_z),\\ \mathbf{H}=|\mathbf{J}|(\mathbf{f}{\zeta }_x+\mathbf{g}{\zeta }_y+\mathbf{h}{\zeta }_z),\end{array}$$

$$\mathbf{J}=\frac{\partial (x,y,z)}{\partial (\xi ,\eta ,\zeta )},\text{and}|\mathbf{J}|=\det (\mathbf{J}).$$

$$\{\begin{array}{l}\tilde{\mathbf{F}}={\mathbf{F}}^l+{\mathbf{F}}^c,\\ \tilde{\mathbf{G}}={\mathbf{G}}^l+{\mathbf{G}}^c,\\ \tilde{\mathbf{H}}={\mathbf{H}}^l+{\mathbf{H}}^c,\end{array}$$

$$\frac{\partial \mathbf{q}}{\partial t}+\frac{\partial {\mathbf{f}}^l}{\partial x}+\frac{\partial {\mathbf{g}}^l}{\partial y}+\frac{\partial {\mathbf{h}}^l}{\partial z}+\frac{1}{|\mathbf{J}|}\left(\frac{\partial {\mathbf{F}}^c}{\partial \xi }+\frac{\partial {\mathbf{G}}^c}{\partial \eta }+\frac{\partial {\mathbf{H}}^c}{\partial \zeta }\right)=0$$

$$\{\begin{array}{l}{\mathbf{F}}^c(\xi ,\eta ,\zeta )=(\tilde{\mathbf{F}}(-1,\eta ,\zeta )-{\mathbf{F}}^l(-1,\eta ,\zeta ))g_L(\xi )+(\tilde{\mathbf{F}}(1,\eta ,\zeta )-{\mathbf{F}}^l(1,\eta ,\zeta ))g_R(\xi ),\\ {\mathbf{G}}^c(\xi ,\eta ,\zeta )=(\tilde{\mathbf{G}}(\xi ,-1,\zeta )-{\mathbf{G}}^l(\xi ,-1,\zeta ))g_L(\eta )+(\tilde{\mathbf{G}}(\xi ,1,\zeta )-{\mathbf{G}}^l(\xi ,1,\zeta ))g_R(\eta ),\\ {\mathbf{H}}^c(\xi ,\eta ,\zeta )=(\tilde{\mathbf{H}}(\xi ,\eta ,-1)-{\mathbf{H}}^l(\xi ,\eta ,-1))g_L(\zeta )+(\tilde{\mathbf{H}}(\xi ,\eta ,1)-{\mathbf{H}}^l(\xi ,\eta ,1))g_R(\zeta ),\end{array}$$

$${\mathbf{f}}_{\mathbf{n},inv}^{com}=\frac{{\mathbf{f}}_{\mathbf{n},inv}^{+}+{\mathbf{f}}_{\mathbf{n},inv}^{-}}{2}-\mathbf{R}|\Lambda |{\mathbf{R}}^{-1}\frac{{\mathbf{q}}^{+}-{\mathbf{q}}^{-}}{2},$$

$$\{\begin{array}{l}{\mathbf{F}}^{num}=|\mathbf{J}||\nabla \xi |{\mathbf{f}}_{\mathbf{n}}^{com}\mathrm{sign}(\mathbf{n}\cdot \nabla \xi ),\\ {\mathbf{G}}^{num}=|\mathbf{J}||\nabla \eta |{\mathbf{f}}_{\mathbf{n}}^{com}\mathrm{sign}(\mathbf{n}\cdot \nabla \eta ),\\ {\mathbf{H}}^{num}=|\mathbf{J}||\nabla \zeta |{\mathbf{f}}_{\mathbf{n}}^{com}\mathrm{sign}(\mathbf{n}\cdot \nabla \zeta ).\end{array}$$

$${\mathbf{q}}^{com}=\frac{{\mathbf{q}}^{+}+{\mathbf{q}}^{-}}{2}.$$

$$\nabla {\mathbf{q}}^{com}=\frac{\nabla {\mathbf{q}}^{+}+{\mathbf{r}}^{+}+\nabla {\mathbf{q}}^{-}+{\mathbf{r}}^{-}}{2},$$

$$\mathbf{r}=\gamma ({\mathbf{q}}^{com}-{\mathbf{q}}_{L/R})\mathbf{n},\text{and}\gamma =|\nabla \varpi |g_{L/R}^{\prime }(\varpi )\mathrm{sign}(\mathbf{n}\cdot \nabla \varpi ),$$

$$\frac{\partial \mathbf{q}}{\partial t}=\mathbf{R}(\mathbf{q}),$$

$${\mathbf{q}}^{n+1}=\Delta t\omega \mathbf{R}({\mathbf{q}}^{n+1})+\sum _{j=1}^s{a}_j{\mathbf{q}}^{n+1-j},$$

$$\{\begin{array}{l}{\mathbf{q}}^{n+1}={\mathbf{q}}^n+\Delta t{\sum }_{i=1}^s{b}_i\mathbf{R}({\mathbf{q}}_i),\\ {\mathbf{q}}_i={\mathbf{q}}^n+\Delta t{\sum }_{j=1}^i{a}_{ij}\mathbf{R}({\mathbf{q}}_j),i=1,\dots ,s,\end{array}$$

$$a_{ii}=\{\begin{array}{ll}0,& i=1,\\ \omega ,& i\ne 1.\end{array}$$

$$\{\begin{array}{l}{\mathbf{q}}^{n+1}={\mathbf{q}}^n+\Delta t{\sum }_{i=1}^s{b}_i\mathbf{R}({\mathbf{q}}_i),\\ {\mathbf{q}}_1={\mathbf{q}}^n,\\ {\mathbf{q}}_i=\Delta t\omega \mathbf{R}({\mathbf{q}}_i)+{\mathbf{q}}^n+\Delta t{\sum }_{j=1}^{i-1}{a}_{ij}\mathbf{R}({\mathbf{q}}_j),i=2,\dots ,s.\end{array}$$

$$\{\begin{array}{l}{\mathbf{q}}^{n+1}={\mathbf{q}}^n+{\sum }_{j=1}^s{m}_j{\mathbf{Y}}_j,\\ {\left(\frac{\mathbf{I}}{\omega \Delta t}-\frac{\partial \mathbf{R}}{\partial \mathbf{q}}\right)}^n{\mathbf{Y}}_i=\mathbf{R}\left({\mathbf{q}}^n+{\sum }_{j=1}^{i-1}{a}_{ij}{\mathbf{Y}}_j\right)+\frac{1}{\Delta t}{\sum }_{j=1}^{i-1}{c}_{ij}{\mathbf{Y}}_j,i=1,2,\dots ,s,\end{array}$$

$$\mathbf{F}({\mathbf{q}}^{\ast })=0,$$

$$\mathbf{F}({\mathbf{q}}^{n+1})=\left(\frac{1}{\omega \Delta t}{\mathbf{q}}^{n+1}-\mathbf{R}({\mathbf{q}}^{n+1})\right)-\frac{1}{\omega \Delta t}\sum _{j=1}^s{a}_j{\mathbf{q}}^{n+1-j}.$$

$$\mathbf{F}({\mathbf{q}}^i)=\left(\frac{1}{\omega \Delta t}{\mathbf{q}}_i-\mathbf{R}({\mathbf{q}}_i)\right)-\frac{1}{\omega \Delta t}\left({\mathbf{q}}^n+\Delta t\sum _{j=1}^{i-1}{a}_{ij}\mathbf{R}({\mathbf{q}}_j)\right),i=2,\dots ,s.$$

$$\frac{{\mathbf{q}}^{k+1,\ast }-{\mathbf{q}}^{k,\ast }}{\Delta \tau }=-\mathbf{F}({\mathbf{q}}^{k+1,\ast }),$$

$${\left(\frac{\mathbf{I}}{\Delta \tau }+\frac{\partial \mathbf{F}}{\partial \mathbf{q}}\right)}^k\Delta {\mathbf{q}}^{k,\ast }=-\mathbf{F}({\mathbf{q}}^{k,\ast }).$$

$${\left(\frac{\partial \mathbf{F}}{\partial \mathbf{q}}\right)}^k={\left(\frac{\mathbf{I}}{\omega \Delta t}-\frac{\partial \mathbf{R}}{\partial \mathbf{q}}\right)}^k.$$

$${\left(\frac{\mathbf{I}}{\Delta \tau }+\frac{\mathbf{I}}{\omega \Delta t}-\frac{\partial \mathbf{R}}{\partial \mathbf{q}}\right)}^k\Delta {\mathbf{q}}^{k,\ast }=-\mathbf{F}({\mathbf{q}}^{k,\ast }).$$

$${\tau }^0={\tau }_{init},\Delta {\tau }^{k+1}=\min \left(\Delta {\tau }^k\frac{||\mathbf{F}|{|}_{L_2}^{k-1}}{||\mathbf{F}|{|}_{L_2}^k},\Delta {\tau }_{max}\right).$$

$$AX=b,$$

$$A=D(\Delta t,\Delta \tau )-\frac{\partial R}{\partial q}.$$

Full text

A comparative study of implicit

Jacobian-free Rosenbrock-Wanner, ESDIRK and BDF methods for unsteady flow simulation with high-order flux reconstruction formulations

Lai Wang  and Meilin Yu PhD candidate. Email: bx58858@umbc.eduAssistant Professor. Corresponding Author. Email: mlyu@umbc.edu Department of Mechanical Engineering,

University of Maryland, Baltimore County, Baltimore, MD 21250

Abstract

We conduct a comparative study of the Jacobian-free linearly implicit Rosenbrock-Wanner (ROW) methods, the explicit first stage, singly diagonally implicit Runge-Kutta (ESDIRK) methods, and the second-order backward differentiation formula (BDF2) for unsteady flow simulation using spatially high-order flux reconstruction/correction procedure via reconstruction (FR/CPR) formulations. The pseudo-transient continuation is employed to solve the nonlinear systems resulting from the temporal discretizations with ESDIRK and BDF2. A Jacobian-free implementation of the restarted generalized minimal residual method (GMRES) solver is employed with a low storage element-Jacobi preconditioner to solve linear systems, including those in linearly implicit ROW methods and those from linearization of the nonlinear systems in ESDIRK and BDF2 methods. We observe that all ROW and ESDIRK schemes (from second order to fourth order) are more computationally efficient than BDF2, and ROW methods can potentially be more efficient than ESDIRK methods. However, the convergence tolerance of the GMRES solver for ROW methods needs to be sufficiently tight to preserve the nominal order of accuracy. In general, ESDIRK methods allow a larger physical time step size for unsteady flow simulation than ROW methods do.

Key Words

Linearly Implicit Rosenbrock-Wanner; ESDIRK; Backward Differentiation Formula; Jacobian-Free; Flux Reconstruction/Correction Procedure via Reconstruction; High-Order Spatiotemporal Methods; Unsteady Flows

1 Introduction

High-order computational fluid dynamics (CFD) methods have been attracting much research attention in past decades due to their superior numerical properties that enable high-fidelity simulation of intricate flows. High-order spatial discretization methods, such as discontinuous Galerkin (DG) methods [1, 2, 3], and flux reconstruction/correction procedure via reconstruction (FR/CPR) methods [4, 5, 6, 7], have shown their capabilities of dealing with turbulent flows [8, 9, 10, 11, 12, 13, 14, 15]. Usually, the high-order explicit strong stability preserving Runge-Kutta (SSPRK) methods [16] are used to integrate the semi-discretized governing equations for unsteady flow simulation. However, due to the Courant–Friedrichs–Lewy (CFL) number constraint, explicit time integration methods may not be the optimal choice for efficient numerical simulation of stiff flow problems, such as wall-bounded turbulent flows at high Reynolds numbers.

The implicit time integration methods have a better stability property than explicit ones. The backward differentiation formula (BDF) is very popular for its ease of implementation. However, BDF methods are not A-stable when the order of accuracy exceeds two. An alternative is the implicit Runge-Kutta (IRK) method [17]. We note that fully coupled IRK methods are not widely used in the CFD community due to the complication of solving fully coupled nonlinear systems. Solution strategies based on the dual time stepping procedure [18] have recently been reported by Jameson [19]. In that work, several fully coupled IRK methods have been evaluated for unsteady flow simulation. To decrease computational complexity, diagonally implicit Runge-Kutta (DIRK) and singly diagonally implicit Runge-Kutta (SDIRK) can be used. A comprehensive review of DIRK methods by Kennedy and Carpenter can be found in Ref. [20]. A recent work by Vermeire and Vincent [21] have analyzed dispersion and dissipation properties of fully-discrete high-order FR methods with two SDIRK schemes, and provided insights on their suitability for implicit large eddy simulation. As a special case of SDIRK, the explicit first stage singly diagonally implicit Runge-Kutta method (ESDIRK) reduces the degree of the nonlinear systems of SDIRK by one. Comparisons of BDF and ESDIRK methods have been performed in Refs. [22, 23]. ESDIRK methods are found to be more efficient than BDF methods. We note that for BDF and DIRK methods, a nonlinear system needs to be solved at each step or each stage.

Rosenbrock methods linearize the nonlinear equations in DIRK methods. Therefore, only one linear equation needs to be solved at each stage. Besides, the Jacobian matrix for this type of multistage methods only needs to be calculated once at the first stage. All these features can potentially make linearly implicit Rosenbrock methods more computationally efficient than DIRK methods. In traditional Rosenbrock methods, exact Jacobian matrices [24, 25] need to be calculated to ensure accuracy. In many stiff flow problems, analytical Jacobian is not easy to obtain; and matrix-based implementation can consume tremendous memory, thus impeding efficient simulation of large-scale flow problems. Rosenbrock-Wanner (ROW) methods  [26, 27, 28] can preserve the nominal order of accuracy with an approximate Jacobian matrix. This flexibility makes it possible to implement a Jacobian-free Krylov subspace solver. The Rosenbrock-Krylov (ROK) methods [29] reformulate the Rosenbrock/ROW method such that stage vectors are obtained from the Krylov subspace using the modified Arnoldi iteration. The ROK methods naturally favor the matrix-free implementation. We notice that even stiffly accurate Rosenbrock methods can suffer from order reduction for moderately stiff problems and improvements have been developed in Ref. [30]. Recent research on the performance of Rosenbrock methods on solving the Navier-Stokes equations with high-order DG-type schemes have been reported in a series of works [31, 32, 33, 34].

If Jacobian-free implementation is used, the advantage that the Jacobian matrix of Rosenbrock methods only needs to be evaluated once for each time step does not hold anymore. In fact, the performance of ROW methods and ESDIRK methods highly depends on how the linear and nonlinear systems are solved. Blom et al. [35] have shown that ROW methods are not necessarily more efficient than ESDIRK methods when a second-order central finite volume method is used. Liu et al. [32] has conducted a comparative study of some third-order ROW methods and a third-order ESDIRK (ESDIRK3) [22] method with a third-order hierarchical WENO (weighted essentially non-oscillatory) reconstructed discontinuous Galerkin (rDG) method. They observed that the third-order ROW methods tested are more efficient than ESDIRK3. Sarshar et al. [36] conducted a comparative study of various matrix-free implicit time-stepping methods on solving two-dimensional (2D) Navier-Stokes equations including SDIRK, Rosenbrock, ROW, and ROK using a second-order spatial discretization. It is found that ROK can be a competitive method compared to other implicit time integrations.

Contributions. In most numerical simulation of the unsteady Navier-Stokes equations with the method of lines, numerical dissipation and dispersion are primarily emphasized for spatial discretization methods, assuming that the influence of time integration methods on numerical simulation is negligible compared to that of spatial discretization. However, the assumption may not hold especially when large time steps, allowed by implicit time integrators, are employed. Several recent works [37, 21, 38] have confirmed that time integration methods can have significant impact on unsteady flow simulation with high-order spatial formulations. This motivates the present study to conduct a systematic comparison of accuracy and efficiency between several widely used implicit time integration methods, including Jacobian-free ROW (second-, third- and fourth-order schemes), ESDIRK (second-, third- and fourth-order schemes) and BDF2, for unsteady numerical simulation of stiff flow problems. In this study, the high-order FR methods are used to discretize the spatial domain. For ESDIRK and BDF2, the pseudo-transient continuation [18] is used to solve the nonlinear systems. The restarted GMRES solver [39] with the element-Jacobi preconditioner serves as the linear solver. The Jacobian-free implementation decreases the tremendous memory consumption of matrix-based implementation. The accuracy and efficiency of the implicit time integrators are compared using several laminar and turbulent flows, including 2D isentropic vortex propagation, 2D laminar flow over a cylinder, and three-dimensional (3D) Taylor-Green vortex evolution.

Article Organization. The rest of the paper is organized as follows. Section 2 gives a brief introduction of the governing equations. The FR method is reviewed in Section 3. Section 4 introduces all the time integration methods tested in this study. Numerical results and corresponding discussions are documented in Section 5. We summarize this work in the last section.

2 Governing Equations

The 3D unsteady compressible Navier-Stokes equations can be written as

[TABLE]

where $\boldsymbol{q}=(\rho,\rho u_{j},E)^{T},\ j=1,2,3,$ is the vector of conserved variables, and $\mathbb{F}$ is the corresponding flux tensor. Specifically, $\rho$ is the fluid density, $u_{j},j=1,2,3,$ are the velocities in three orthogonal directions, and $E=\frac{p}{\gamma-1}+\frac{1}{2}\rho\sum_{k=1}^{3}(u_{k}u_{k})$ is the total energy per unit volume, where $p=\rho RT$ is the pressure, $T$ is the temperature, $R=C_{p}-C_{v}$ is the ideal gas constant, $\gamma$ is the specific heat ratio defined as $\gamma=C_{p}/C_{v}$, and $C_{p}$ and $C_{v}$ are specific heat capacities at constant pressure and volume, respectively. In this study, $\gamma$ is set as 1.4. The flux tensor $\mathbb{F}$ consists of the inviscid part and viscous part, which can be expressed as $\mathbb{F}=\mathbb{F}_{inv}(\boldsymbol{q})-\mathbb{F}_{vis}(\boldsymbol{q},\nabla\boldsymbol{q})$. Note that $\mathbb{F}_{inv}(\boldsymbol{q})$ and $\mathbb{F}_{vis}(\boldsymbol{q},\nabla\boldsymbol{q})$ can be rewritten in a vector format as: for any $j,\ j\in\{1,2,3\},$

[TABLE]

where $\delta_{ij}$ is the Kronecker delta, $\tau_{ij}$ is the viscous stress defined as

[TABLE]

and following Fourier’s law, the heat flux $K_{j}$ is defined as $K_{j}=-\kappa\partial T/\partial x_{j}$. In this study, the thermal conductivity $\kappa$ is calculated from the Prandtl number $Pr$ as $\kappa=\mu C_{p}/Pr$. The fluid viscosity $\mu$ and the Prandtl number $Pr$ are treated as constants.

3 Flux Reconstruction Methods

To achieve an efficient implementation of Eq. (1), we transfer the Navier-Stokes equations from the physical domain $(x,y,z)$ to the computational domain $(\xi,\eta,\zeta)$. Thus Eq. (1) can be expressed as

[TABLE]

where

[TABLE]

and

[TABLE]

Herein, $\boldsymbol{f}=\mathbb{F}_{1}$, $\boldsymbol{g}=\mathbb{F}_{2}$, and $\boldsymbol{h}=\mathbb{F}_{3}$ are flux vectors in the $x$, $y$, and $z$ directions of a Cartesian coordinate system. ${\partial\xi_{i}}/{\partial x_{j}},\ i,j=1,2,3,$ are metrics in the coordinate transformation.

The flux polynomial reconstructed from the FR method consists of two parts, one of which is the local flux polynomial and the other is the correction polynomial. On solving Eq. (3), the reconstructed polynomials $\widetilde{\boldsymbol{F}}$, $\widetilde{\boldsymbol{G}}$ and $\widetilde{\boldsymbol{H}}$ of $\boldsymbol{F}$, $\boldsymbol{G}$ and $\boldsymbol{H}$ can be expressed as

[TABLE]

where the superscript ‘$l$’ stands for the local flux and ‘$c$’ stands for the correction flux. Consequently, Eq. (1) can be rewritten as

[TABLE]

For hexahedral elements, $\boldsymbol{F}^{c}$ , $\boldsymbol{G}^{c}$ and $\boldsymbol{H}^{c}$ can be explicitly expressed as

[TABLE]

where $g_{L/R}$ are the correction polynomials. In this study, we employ the Radau polynomials to recover the nodal FR-DG method. $\widetilde{\boldsymbol{F}}$, $\widetilde{\boldsymbol{G}}$ and $\widetilde{\boldsymbol{H}}$ at element interfaces are referred as numerical fluxes $\boldsymbol{F}^{num}$, $\boldsymbol{G}^{num}$ and $\boldsymbol{H}^{num}$. The inviscid common fluxes can be obtained from approximate Riemann solvers. In this study, the Roe approximate Riemann solver [40] is used to calculate the common fluxes at the cell interfaces in their normal directions as

[TABLE]

where superscripts ‘$-$’ and ‘$+$’ denote the left of right side of the current interface, the subscript $\boldsymbol{n}$ is the unit normal direction from left to right, $\boldsymbol{f}_{\boldsymbol{n}}=fn_{x}+gn_{y}+hn_{z}$ is the flux projection in the $\boldsymbol{n}$ direction, $\Lambda$ is a diagonal matrix consisting of the eigenvalues of the Jacobian $\partial\boldsymbol{f}_{\boldsymbol{n}}/\partial\boldsymbol{q}$, and $\boldsymbol{R}$ consists of the corresponding right eigenvectors evaluated with the Roe-averaged variables. Numerical common fluxes can be obtained as

[TABLE]

The common viscous fluxes at the cell interfaces are $\boldsymbol{f}_{n,vis}^{com}=\boldsymbol{f}_{vis}(\boldsymbol{q}^{+},\nabla\boldsymbol{q}^{+},\boldsymbol{q}^{-},\nabla\boldsymbol{q}^{-})$. Here we need to define the common solution $\boldsymbol{q}^{com}$ and common gradient $\nabla\boldsymbol{q}^{com}$ at the cell interface. By simply taking average of the primitive variables, we get

[TABLE]

The common gradient is computed as

[TABLE]

where $\boldsymbol{r}^{+}$ and $\boldsymbol{r}^{-}$ are the corrections to the gradients on the interface. The second approach of Bassi and Rebay (BR2) [8] is used to calculate the corrections. For the hexahedral element, the correction terms are defined as

[TABLE]

where $\varpi\in\{\xi,\eta,\zeta\}$, and $\boldsymbol{q}_{L/R}$ is the local solution on the interface. If the interface is the left boundary of the element, then the local solution $\boldsymbol{q}_{L}$ is used, and $\varpi$ is $-1$; if the interface is the right boundary of the element, then the local solution $\boldsymbol{q}_{R}$ is used, and $\varpi$ is $1$. In this study, $g(\pm 1)=\pm(P+1)(P+2)/2$ is used to stabilize the FR method  [41], where $P$ is the polynomial degree of the solution.

4 Time Integration Methods

To solve Eq. (7), we rewrite this equation in a semi-discretized form as

[TABLE]

for the convenience of notation. We use the method of lines approach to solve this equation: the partial differential equation is first discretized in space with the FR formulation in Section 3, and then is marched in time with the time integrators described below.

4.1 The backward differentiation formula

The BDF method for solving Eq. (14) can be expressed as

[TABLE]

BDF methods are not A-stable when the order of accuracy exceeds two. In this study, we only consider the BDF2 method for comparison. For BDF2, $s=2$, $\omega=2/3$, $a_{1}=4/3$ and $a_{2}=-1/3$.

4.2 The explicit first stage, single diagonally implicit Runge-Kutta method

The ESDIRK methods can be written as

[TABLE]

where $s$ is the number of stages. The second-order, three-stage ESDIRK2 [20], third-order, four-stage ESDIRK3 [22] and fourth-order, six-stage ESDIRK4 [22] methods are studied in this paper. All the ESDIRK methods investigated in this study have the feature that

[TABLE]

Therefore, we can rewrite Eq. (16) as

[TABLE]

4.3 Linearly implicit Rosenbrock methods

The general form of Rosenbrock methods for solving Eq. (14) from time step $n$ to $n+1$ can be written as

[TABLE]

where $s$ is the number of stages. Rosenbrock-Wanner methods have the flexibility of evaluating the Jacobian matrix approximately while preserving the accuracy of the schemes. In this study, we restrict our attention to three popular stiffly accurate Rosenbrock-Wanner methods, i.e., the second-order, three-stage ROS2PR [27], the third-order, four-stage ROS34PW2 [26] and the fourth-order, six-stage RODASP [28]. We note that it has been shown in Ref. [30] that ROS34PW2 and RODASP can suffer from order reduction. To simplify the notation, we use ROW2, ROW3 and ROW4 to denote them, respectively.

4.4 Iterative methods

The nonlinear equations in Eq. (15) and Eq. (18) can be written as

[TABLE]

where $\boldsymbol{q}^{*}$ is $\boldsymbol{q}^{n+1}$ for BDF and $\boldsymbol{q}^{i}$ for ESDIRK, respectively. For BDF methods,

[TABLE]

For ESDIRK methods,

[TABLE]

In this work, the pseudo-transient continuation is employed to solve the nonlinear equations [18, 22, 19]. The pseudo-transient continuation is an alternative of inexact Newton’s method to solve the steady state equation $\boldsymbol{F}(\boldsymbol{q}^{*})=0$ iteratively as

[TABLE]

Herein, $k$ is the iteration step for the pseudo-transient continuation. Eq. (23) can be linearized as

[TABLE]

For a steady problem, as $k\to\infty$, $\tau\to\infty$ and $\Delta\boldsymbol{q}^{k}\to 0$. Therefore, $\boldsymbol{F}(\boldsymbol{q}^{k,*})\to\boldsymbol{F}(\boldsymbol{q}^{*})$. The Jacobian matrix for fully implicit methods can be expressed as

[TABLE]

Substitute the Jacobian matrix into Eq. (24) to obtain the final form of the linear system in the pseudo-transient continuation procedure as follows

[TABLE]

The pseudo-transient continuation procedure requires an adaptation algorithm of the pseudo time step size to complete the method. In this study, we employ successive evolution relaxation (SER) algorithm [42] as

[TABLE]

As the pseudo time marches forward, a series of linear equations (26) are successively solved until convergence. Ideally, we would expect that when $\Delta\tau$ approaches $\Delta\tau_{max}=\infty$, the residual of the pseudo-transient procedure will gradually converge to machine zero. However, in our practice of simulating wall-bounded flows with grids clustered in near wall regions, we have to choose a moderately large $\Delta\tau_{max}$ to ensure that the residual of the linear solver (e.g., the restarted GMRES as will be explained in the next subsection) can at least drop by one magnitude when $\Delta\tau$ equals to $\Delta\tau_{max}$. Once the linear solver fails, we would reject the current pseudo-transient continuation iteration and decrease $\Delta\tau_{max}$ by half to continue. In this study, if not specifically mentioned, $\Delta\tau_{init}=\Delta t$ and $\Delta\tau_{max}=10^{20}$.

4.5 Jacobian-free implementation of GMRES

The restarted GMRES is employed to solve Eq. (26) as well as the linear system in Eq.(19). All linear systems can be expressed as

[TABLE]

where

[TABLE]

The first term $D(\Delta t,\Delta\tau)$ on the right-hand side of the above equation is a diagonal matrix related to $\Delta t$ and $\Delta\tau$. The GMRES method approximates the exact solution by a vector $x_{n}\in K_{n}$ that minimizes the Euclidean norm $||Ax_{n}-b||$ where $K_{n}$ is the $n$-th Krylov subspace

[TABLE]

In the GMRES solver, the sparse matrix $\frac{\partial\boldsymbol{R}}{\partial\boldsymbol{q}}$ in $A$ only appears in the matrix-vector product. For an unknown vector $\boldsymbol{X}$, the matrix-vector product can be approximated as

[TABLE]

$\epsilon=10^{-8}$ in this study. Interested readers are referred to Ref. [43] for more discussions on this approximation.

The performance of the Newton-Krylov methods substantially depends on the preconditioner. In the context of Jacobian-free implementation, the diagonal blocks of the Jacobian matrix, i.e., the element-Jacobi preconditioner, can be an effective choice among low-storage preconditioners, such as the matrix-free LU-SGS (Lower-Upper Symmetric-Gauss-Seidel) preconditioner [44], and $p$-multigrid preconditioner [45]. In this study, only the element-Jacobi preconditioner is considered. The restart number of the GMRES solver is set as 60.

5 Numerical Results

All the simulations presented in this section are performed on the High Performance Computing Facility (HPCF) of the University of Maryland, Baltimore county (UMBC). All nodes used in simulations have two Intel E5-2650v2 Ivy Bridge (2.6 GHz, 20 MB cache) processors with eight cores apiece, for a total of 16 cores per node. A quad data rate (QDR) Infiniband switch connects all the nodes. Ideally the system can achieve a latency of 1.2 $\mu$sec to transfer a message between two nodes, and can support a bandwidth of up to 40 gigabits per second (40Gbps). Every node possesses 64 GB RAM. All nodes are running Red Hat Enterprise Linux 6.4. We employ g++ (GCC) 4.8.4 with mpich-3.1.4 to compile the code for parallel simulation.

5.1 Vortex propagation

The 2D isentropic vortex propagation problem is employed as a benchmark test test case to investigate the accuracy and efficiency of different time integrators in this subsection. The free stream flow conditions are set as $(\rho,u,v,Ma)=(1,1,1,0.5)$ and the ideal gas constant $R$ is set as $1.0$ for this case. The perturbation is defined as [31]

[TABLE]

where $\phi=\frac{1}{2}$ and $\alpha=5$ are parameters that define the vortex strength, and $r=(x-x_{0})^{2}+(y-y_{0})^{2}$ is the distance to the center of the vortex $(x_{0},y_{0})=(0,0)$ at $t=0$. The periodic domain is defined as $\Omega=[-10,10]^{2}$. We simulate this problem on a uniform mesh of $50\times 50$ elements with the $P^{6}$ (i.e. $7th$ order) FR method for one period. The time step size $\Delta t$ is refined from $T/100$ to $T/800$. The error of any variable $s$ is defined as

[TABLE]

where $s_{exact}$ is the exact value, $s_{num}$ is the numerical value, and $V$ is the volume of the domain $\Omega$.

To monitor the convergence of the nonlinear system (20) and linear system (28), two convergence criteria are set up: one is for pseudo time iterations, denoted as $tol_{rel,nonlinear}$; the other is for the restarted GMRES linear solver, denoted as $tol_{rel,linear}$. Herein, the subscript ‘$rel$’ stands for ‘relative’. This is due to that the relative residual with respect to that at the first step (i.e., the residual at the first step in pseudo time iterations for a nonlinear system, and the residual at the first step in the restarted GMRES for a linear system) is used in this study.

5.1.1 Effect of pseudo-transient continuation convergence criterion and GMRES convergence criterion on ESDIRK and BDF2

The pseudo-transient continuation does not require exact solution of the linear system at each iteration. We first employ a relative big tolerance $tol_{rel,linear}=10^{-1}$ for GMRES to investigate the effect of $tol_{rel,nonlinear}$ on accuracy and efficiency of fully implicit methods. The convergence study of ESDIRK and BDF2 of different $tol_{rel,nonlinear}$, i.e., $10^{-2}$, $10^{-4}$, $10^{-6}$ and $10^{-8}$, is presented in Table LABEL:fullyimplicit_nonlinear. We observe that when $tol_{rel,nonlinear}$ is refined from $10^{-2}$ to $10^{-4}$, all ESDIRK methods studied here will preserve the nominal order of accuracy except that ESDIRK4 shows slightly order reduction when $\Delta t$ is refined from $T/400$ to $T/800$. ESDIRK methods have shown better accuracy and efficiency than BDF2. Generally, the higher the order of accuracy of ESDIRK is, the more efficient ESDIRK is (see Figure 1).

We have also conducted a study on the effect of $tol_{rel,linear}$ on efficiency when a small enough $tol_{rel,nonlinear}=10^{-6}$ is employed. $tol_{rel,linear}$ spans in $\{10^{-1},10^{-2},10^{-4},10^{-6}\}$. The timing results are presented in Table 2. For this specific problem, all values of $tol_{rel,linear}$ lead to the same numerical errors as expected since the accuracy of ESDIRK and BDF2 methods solely depends on the convergence of pseudo-transient continuation. As documented in Table 2, the computational cost will keep on increasing when we refine the tolerance $tol_{rel,linear}$. This indicates that in order to save computational cost, $tol_{rel,linear}$ for ESDIRK and BDF2 needs not to be tight.

5.1.2 Effect of GMRES convergence criterion on ROW

A comparison of different convergence criteria of the restarted GMRES solver is conducted to study its impact on the accuracy and efficiency of the ROW time integrators. It is observed that when the convergence criterion is not tight, such as $tol_{rel,linear}=10^{-2}$ and $10^{-4}$, the ROW methods cannot preserve the nominal order of accuracy. This is not a surprise considering that the residual convergence is directly related to the accuracy of the solution in ROW. When $tol_{rel,linear}$ is tight, such as $tol_{rel,linear}=10^{-6}$ and $10^{-8}$, all ROW methods can preserve the nominal order of accuracy, excepted that ROW4 shows order reduction when the time step is refined from $T/400$ to $T/800$. We also notice that when $tol_{rel,linear}$ is refined from $10^{-6}$ to $10^{-8}$, no significant differences in errors are observed. Another observation is that when $tol_{rel,linear}$ is overrefined, the CPU time of the simulation is noticeably increased. In general, we do not recommend machine zero convergence criterion for $tol_{rel,linear}$. A relatively tight value such as $10^{-6}$ is sufficient to preserve the accuracy of high-order Rosenbrock methods for this problem.

5.1.3 Comparison of different time integrators

The convergence study of ESDIRK and BDF2 and linearly implicit ROW are summarized in Figure 1. For ESDIRK and BDF2, $tol_{rel,nonlinear}^{ESDIRK,BDF2}=10^{-6}$ and $tol_{rel,linear}^{ESDIRK,BDF2}=10^{-1}$. For Rosenbrock methods, $tol_{rel,linear}^{ROW}=10^{-6}$.

Figure 1(b) presents the CPU time versus errors of different time integrators with different convergence criteria. As illustrated in Figure 1(b), all multistage methods are significantly more efficient than BDF2. We notice that ROW2 and ESDIRK2 intersect with ROW3 and ESDIRK3. However, as the error threshold is decreased, higher-order methods will be more efficient.

Generally speaking, when $tol_{rel,nonlinear}^{ESDIRK}$ is the same as $tol_{rel,linear}^{ROW}$, ROW methods are more efficient than ESDIRK methods. However, ROW methods cannot preserve the nominal order of accuracy when the convergence criterion is not tight, and even suffer from severe order reduction [30]. Instead, $tol_{rel,nonlinear}=10^{-4}$ can make ESDIRK methods preserve the nominal order of accuracy (0.46 order reduction at most). It is observed that when $tol_{rel,linear}^{ROW}=10^{-6}$ and $tol_{rel,nonlinear}^{ESDIRK}=10^{-4}$ are employed, ESDIRK methods are more efficient than ROW methods. This indicates that ESDIRK methods tend to be over-solved more easily than ROW methods if the nonlinear convergence criterion $tol_{rel,nonlinear}^{ESDIRK}$ is not set up judiciously. We have also noticed that ESDIRK can be more robust than ROW. As documented in Table LABEL:fullyimplicit_nonlinear and Table LABEL:linearly_implicit, when the tolerance criteria are set to $10^{-2}$, ESDIRK2 can show optimal convergence rate, while ROW cannot preserve it; when $\Delta t=T/100$, ROW2 even fails to converge.

5.2 Laminar flow over a circular cylinder

In this section, we employ the laminar flow over the circular cylinder as an example to study the performance of different time integrators. This case has been tested in various literature [22, 35]. The flow conditions are set as $Ma=0.1$ and $Re=1200$, where $Ma$ stands for Mach number, and $Re$ stands for Reynolds number. The diameter of the cylinder is set as $1$, and computational domain is $[-100,200]\times[-100,100]$. The mesh in the near wall region and the wake region, and an instance of the vortex shedding are presented in Figure 2. There are 5690 elements in the mesh. The height of the first layer of the mesh is roughly $0.0033$.

The $P^{3}$ FR method is employed for the spatial discretization. $\Delta\tau_{init}=0.01$ and $\Delta\tau_{max}={1.0}$ are used for all simulations in this section. A moderate $\Delta\tau_{max}$ is to prevent that GMRES solver can not even converge by one magnitude when $\Delta\tau=\Delta\tau_{max}$. $tol_{rel,linear}^{ROW}=tol_{rel,nonlinear}^{ESDIRK,BDF2}=10^{-6}$, and $tol_{rel,linear}^{ESDIRK,BDF2}=10^{-1}$. As aforementioned, for ESDIRK and BDF, if the linear solver fails to drive the residual to drop by one magnitude, the current iteration in the pseudo-transient continuation procedure will decrease $\Delta\tau_{max}$ by half and restart.

The flow is initialized with the steady solution at $Re=40$. And then we use ESDIRK4 to run this simulation untill $t=180$ with $\Delta t=0.001$ to obtain the initial conditions for the convergence and efficiency study. For the convergence and efficiency study, we run all simulations for ten seconds along the physical time. The time step size is refined from $0.2$ to $0.00625$. We use the numerical results of explicit SSPRK3 with a small time step $\Delta t=5\times 10^{-6}$ as the reference value. The drag coefficient $C_{d}$ is used for the error estimation. The error is calculated as

[TABLE]

where $C_{d,ref,n}$ is the reference value from SSPRK3, and $N$ is the number of time steps. The results from convergence and efficiency study are presented in Figure 3(a) and Figure 3(b), respectively. As the time step size is refined, all second- and third-order methods will converge at the nominal convergence rate. For both ROW4 and ESDIRK4, we have observed order reduction. The order reduction of ROW4 is more severe than ESDIRK4. In terms of CPU time, when $\Delta t=0.2$, unexpected computational cost is observed for ESDIRK and BDF methods. Many iterations in the pseudo-transient continuation are rejected due to the poor performance of the element-Jacobi preconditioner. However, all ROW methods fail when $\Delta t=0.2$ and ROW2 even fails when $\Delta t=0.1$. At a relatively large error level, such as $10^{-3}$, second-order methods take the least mount of time. However, to reach a lower error level, higher-order methods are more efficient. We have noticed that the abnormal increase in CPU time for ROW4 when $\Delta t=0.05$ and $\Delta t=0.025$. This is due to the fact that the residual of the restarted GMRES solver with an element-Jacobi preconditioner sometimes cannot converge to the designated $tol_{rel,linear}$ when the maximum number of iteration is reached. From this study, we find that the performance of schemes from the ESDIRK family is more consistent than that of ROW methods when simulating unsteady flows over walls; and when the requirement of the error threshold is stringent, such as $10^{-6}$ in this case, ROW can be more efficient than ESDIRK.

5.3 Taylor-Green vortex

The Taylor-Green vortex is a benchmark to test the accuracy and performance of high-order methods on the direct numerical simulation of a 3D periodic and transitional flow defined by initial conditions [32]

[TABLE]

The domain is $\Omega=[-\pi L,\pi L]^{3}$. The Reynolds number of the flow is defined as $Re=\frac{\rho_{0}V_{0}L}{\mu}$ and is equal to 1600. For this study, we consider the flow with weak compressibility and the perfect gas law holds, i.e., $p=\rho RT$. The Prandtl number is $Pr=\frac{\mu C_{p}}{\kappa}=0.71$. We assume that the gas has zero bulk viscosity $\mu_{v}=0$. The Mach number $Ma=\frac{V_{0}}{c_{0}}=0.1$, where $c_{0}$ is the speed of sound corresponding to $p_{0}$. The characteristic convection time is defined as $t_{c}=\frac{L}{V_{0}}$. The maximum dissipation occurs at $t\approx 8t_{c}$. An uniform $64^{3}$ mesh is employed and the $P^{3}$ FR methods is used for simulation. We conduct all simulation until $t=10t_{c}$. Figure 4 presents the isosurface of the Q criterion with $Q_{criterion}=1$ colored by $Ma$ at $t=8t_{c}$.

We employ the error of kinetic energy dissipation rate for the accuracy and efficiency study. The kinetic energy dissipation rate of compressible flows is the summation of three contributions as $\epsilon=\epsilon_{1}+\epsilon_{2}+\epsilon_{3}$:

[TABLE]

where $\boldsymbol{S}^{d}$ is the deviatoric part of the strain rate tensor, and $V$ is the volume of the domain $\Omega$,

[TABLE]

where $\mu_{v}=0$, and

[TABLE]

The numerical results of SSPRK3 method with $\Delta t=2\times 10^{-4}t_{c}$ is adopted as the reference data for error evaluations. We define the error of the kinetic energy dissipation rate as

[TABLE]

where $\epsilon_{ref,n}$ is the reference value from SSPRK3, and $N$ is the number of time steps.

As discovered in previous sections, the convergence criteria have significant effect on the efficiency of ROW, ESDIRK and BDF2. In this section, we only consider $tol_{rel,linear}$ of ROW is the same as $tol_{rel,nolinear}$ of ESDIRK and BDF2. Herein, $tol_{rel,linear}^{ROW}=tol_{rel,nonlinear}^{ESDIRK,BDF2}=10^{-6}$. For the inexactly linear-solving part of ESDIRK and BDF2, we employ $tol_{rel,linear}^{ESDIRK,BDF2}=10^{-1}$ to save computational cost.

The time step size $\Delta t$ is refined from $t_{c}/25$ to $t_{c}/100$. The kinetic energy dissipation rate history when $\Delta t=t_{c}/25$ is presented in Figure 5(a). A close-up view within $t/t_{c}\in[8.5,10]$ is illustrated in Figure 5(b). The numerical results from the spectral method on a $512^{3}$ mesh is also presented for reference [46]. Our observation is that the results of ROW4 and ESDIRK4 almost coincide with that of SSPRK3. The result of BDF2 is the least accurate. The convergence study is presented in Figure 6. Figure 6(a) shows the error vs. time step size $\Delta t/t_{c}$ and Figure 6(b) shows the error vs. CPU time. From Figure 6(a), the convergence features of ROW and ESDIRK are almost the same. As is shown in Figure 6(b), BDF2 is not efficient as expected. When $tol_{rel,linear}^{ROW}=tol_{rel,nonlinear}^{ESDIRK,BDF2}=10^{-6}$, ROW methods are found to be more efficient than ESDIRK methods. From Figure 5, we find that for turbulent simulation, compared to the results from SSPRK3 with very small time steps, all time integrators with excessively large time step size will lead to numerical dissipation of the kinetic energy dissipation rate except ROW4 and ESDIRK4. This reveals that the dissipation due to the temporal disretization should also be taken into account for turbulent simulation. We refers interested readers to Refs. [21, 38] for more discussions.

6 Conclusions

In this study, we compare the accuracy and efficiency of ROW and ESDIRK (from second order to fourth order) and BDF2 for unsteady flow simulation with high-order FR formulations. We find that the multistage time integrators, i.e. ROW and ESDIRK, are more computationally efficient than BDF2. When we compare ROW and ESDIRK, the efficiency of ROW and ESDIRK highly depends on the convergence criteria of solving nonlinear equations and linear equations. For most of the case, when $tol_{rel,linear}^{ROW}=tol_{rel,nonlinear}^{ESDIRK}$, ROW is more efficient. However, since the residual convergence is directly related to the accuracy of the solution with ROW, the tolerance $tol_{rel,linear}^{ROW}$ needs to be tight for unsteady flow simulation. By contrast, the tolerance $tol_{rel,nonlinear}^{ESDIRK}$ for pseudo-time iterations in ESDIRK needs not to be tight, and the tolerance $tol_{rel,linear}^{ESDIRK}$ for the restarted GMRES linear solver can always be set to a large value, i.e. $10^{-1}$ in all cases tested here. Therefore, when $tol_{rel,nonlinear}^{ESDIRK}$ is allowed to be larger than $tol_{rel,linear}^{ROW}$, ESDIRK can be more efficient than ROW. Besides, ROW methods are easier to show order reduction than ESDIRK methods, making their performance less consistent than that of ESDIRK methods when flow problems become complicated, such as unsteady flows over walls. We also find that the performance of the preconditioner can substantially affect the efficiency of ROW methods. Since the element-Jacobi preconditioner is used in this study, time steps for ROW cannot be very large in order to maintain good performance. Instead, ESDIRK methods generally allows a larger time step size for unsteady simulations than ROW methods.

Acknowledgment

The authors gratefully acknowledge the support of the Office of Naval Research through the award N00014-16-1-2735, and the faculty startup support from the department of mechanical engineering at the University of Maryland, Baltimore County (UMBC). The hardware used in the computational studies is part of the UMBC High Performance Computing Facility (HPCF). The facility is supported by the U.S. National Science Foundation through the MRI program (grant nos. CNS-0821258, CNS-1228778, and OAC-1726023) and the SCREMS program (grant no. DMS-0821311), with additional substantial support from UMBC.