Stability evaluation of high-order splitting method for incompressible flow based on discontinuous velocity and continuous pressure

In this work, we deal with high-order solver for incompressible flow based on velocity correction scheme with discontinuous Galerkin discretized velocity and standard continuous approximated pressure. Recently, small time step instabilities have been reported for pure discontinuous Galerkin method, in which both velocity and pressure are discretized by discontinuous Galerkin. It is interesting to examine these instabilities in the context of mixed discontinuous Galerkin–continuous Galerkin method. By means of numerical investigation, we find that the discontinuous Galerkin–continuous Galerkin method shows great stability at the same configuration. The consistent velocity divergence discretization scheme helps to achieve more accurate results at small time step size. Since the equal order discontinuous Galerkin–continuous Galerkin method does not satisfy inf-sup stability requirement, the instability for high Reynolds number flow is investigated. We numerically demonstrate that fine mesh resolution and high polynomial order are required to obtain a robust system. With these conclusions, discontinuous Galerkin–continuous Galerkin method is able to achieve high-order spatial convergence rate and accurately simulate high Reynolds flow. The solver is tested through a series of classical benchmark problems, and efficiency improvement is proved against pure discontinuous Galerkin scheme.


Introduction
Discontinuous Galerkin method (DGM) is one of the most potential high-order discretization method among the state-of-the-art methods, such as finite difference methods and finite volume method (FVM). DGM has attracted lots of attention from both academic and industry community for the easiness to achieve highorder spacial convergence rate, geometrical flexibility, numerical stability, and good extensibility. On the contrary, the numerical simulation of incompressible Navier-Stokes (INS) equations is a key issue in the research of DGM but still far from completeness. This article is devoted to discussing the stability DGM coupled with continuous Galerkin (CG) in the simulation of INS problems.
For the simulation of INS equations, coupled solver solves all components of velocity and pressure in a monolithic manner, which require the solution of saddle point problem and non-linear iteration. Due to the complexity, the coupled solver has only been applied in small scale academic cases. An alternative choice is splitting method proposed by Guermond et al. 1 and Shahbazi et al. 2 to solve the strong coupling of pressure and velocity. Contrast to the coupled solver, a Poisson and Helmholtz equation are solved sequentially for pressure and velocity, which shows to be much more effective in many applications. Splitting methods are also widely used in the context of DGM discretization.
In recent years, several kinds of instabilities have been reported for splitting methods when employed with DGM. Ferrer and Willden 3 show the small time step instability in the simulation of unsteady Stokes and vortex flow with this scheme. They find that a minimum time step Dt(h, k) is required to achieve long-time stable simulation. This idea is extended by Ferrer et al. 4 showing that the minimum time step size is constrained by C h=k h 2 k À3 for both continuous and DGM. A series of papers are published trying to find out the reasons behind and remedy these problems. Krank et al. 5 attribute the instability to the non-zero velocity divergence and a div-div stabilization approach is proposed. Similar work has also been done by Joshi et al. 6 to wipe the spurious divergence at element interfaces. Emamy et al. 7 propose a modified version of projection scheme to reduce the influence of time step size. Steinmoeller et al. 8 also review the instabilities of DG in simulating under-resolved and small viscosity cases. And they suggest the instabilities are related to the inf-sup instability of DGM. Recently, Fehn et al. 9 propose that the small time instability is due to the discretization of velocity divergence and pressure gradient. In the numerical investigation, they find that a stable and robust method can be obtained using integration by parts.
In the theory of finite element, a mixed-order continuous velocity and discontinuous pressure are usually choosen to solve the INS problem, which exhibit good numerical behavior. But considering moderate-to-high Reynolds flow, discontinuous velocity helps to keep stability when convection term becomes dominant. As a supplement to pure DG (both of velocity and pressure are discretized with DGM), Botti and Di Pietro 10 present a pressure-correction scheme for the INS equations, which discrete velocity with standard DGM while pressure with CG. A major preliminary of this work is that the procedure of solving the pressure during the splitting method is the most expensive stage, while approximation of the pressure with CG can significantly reduce computational cost compared to monolithic pure DG. This is much more favorable for three-dimensional (3D) simulations. Meanwhile, the stability for convection dominant problem is also inherited. This scheme has been validated against a large set of classical benchmarks and shows great computational efficiency improvement.
However, the instabilities encountered in pure DG method for small time step size and under-resolved cases may also exhibit in mixed DG-CG method. To our best knowledge, the evaluation on the stability of DG-CG method for INS equations has not been published. In this article, we develop a high-order INS solver based on velocity correction scheme, in which mixed equal order DG-CG method is used. The stabilities of DG-CG method for small time step size and small viscosity are numerically investigated. This work is crucial to the employment of DG-CG method and may also shed light on the remedy of instabilities in pure DG scheme. Major innovations of this article include the following: The behavior of DG-CG method for small time step instability is investigated. On the simulation of Pearson vortex problem, DG-CG method shows great stability on the regimes where pure DG method by Hesthaven and Warburton 11 reported to be unstable. In addition, the consistent velocity divergence discretization scheme helps to achieve more accurate results at small time step size. The spatial stability of DG-CG method in low viscosity situations is numerically investigated through comprehensive experiments for steady state Navier-Stokes equations. We find that sufficient fine mesh and high polynomial order help to make DG-CG scheme more robust and accurate for high Reynolds number flow. The DG-CG solver is validated through a series of classical problems. With the guidance of our conclusions, stabilities for relative high Reynolds flow and efficiency against pure DG scheme have been proved.
The remaining parts of this article are organized as follows: We start by introducing the governing equations and velocity correction steps. In the ''Spatial discretization'' section, the spatial discretization is presented in detail. ''Numerical results'' section numerically investigate the small time and spatial stability through a serious of classical benchmark problems. The temporal and spatial convergence rate and computational efficiency are also examined. Finally, this article in concluded in the ''Conclusion'' section.

Governing equations
Here, we consider the unsteady INS equations in conservation form defined on physical domain O ∂u ∂t + r Á F (u) = À rp + nr 2 u ð1Þ where u = (u 1 , :::, u d ) T is the velocity, p indicates the kinematic pressure, and n is the kinematic viscosity. The definition of convection flux F reads F (u) = u u.
For a particular problem, the initial conditions are given as u(t = 0) = u 0 and p(t = 0) = p 0 in O, where u 0 and p 0 also have to fulfill the governing equations constraints. The Dirichlet conditions defined on ∂O D and Neumann conditions on ∂O N are prescribed as with the constraints ∂O D [ ∂O N = ∂O and ∂O D \ ∂O N = [, respectively. Here,n denotes the outward facial normal vectors. Since the gradient instead of absolute value of pressure makes sense in INS, the resultant value of p can be adjusted by imposing p = g p on ∂O N . In case of pure Dirichlet boundary conditions, two additional constraints are needed, with Ð ∂O g Ándx = 0 to satisfy the continuity equation and Ð O pdx = 0 to guarantee unique pressure solution.

Velocity correction method and temporal discretization
The velocity correction method combined with discontinuous Galerkin (DG) has been investigated by several works, including Hesthaven and Warburton, 11 Ferrer et al., 4 and Emamy et al. 7 Here, the fixed semi-explicit dual-splitting schemes are used to solve equations (1) and (2). Within this scheme, the INS equations are split into three distinct substeps which can be solved successively.
In the first substep, equation (1) is advanced using Adams-Bashforth time scheme with the nonlinear convection term treated explicitly where J = 1, 2 corresponding to different temporal order. The value of constant time integration coefficients can refer to Karniadakis et al. 12 However, this temporal scheme should be performed from J = 1 and increase the order successively, for it can not start automatically when J . 1.
The second substep is pressure projection, in which the intermediate velocity is updated with the following equation By projecting the velocity into weakly divergence free space, equation (6) turns into a Poisson equation for pressure at time t n + 1 as Specially, Neumann boundary conditions should be applied on O D for pressure to keep consistent with governing equations. Multiplying the momentum equation (1) with face normal, we obtain the formulation of pressure gradient given aŝ where v = r 3 u is the vorticity. It should be noted that only the solenoidal part of viscosity is accounted for instead of original formulation, yielding the constraint of divergence free. Then, intermediate velocity can be updated with p n + 1 bỹ The final substep takes the viscosity into consideration and solves an implicit Helmholtz equation Spatial discretization

DG and CG space
The space discretization deployed here is based on the combination of DG and standard CG method. This choice benefits from the stability of DG on handling convection term and efficiency of CG in solving the projection step. The computational domain denoted by O h = [ K k = 1 D k is a numerical approximation of O, which consists of K geometric compatible nonoverlapping cells D k . Theoretically, the shape of D k can be arbitrary. Every interior face ∂D À k \ ∂D + k contains information from two adjacent cells, where the current cell is marked by a symbol ''-'' while adjacent cell by '' + .'' With these notations, the average of a scalar or vector variable u on a face is defined as The discontinuity of variables across adjacent cell faces may be expressed in similar formulation To discrete in either continuous or discontinuous formulation, approximate solution u h in cell k is represented by a linear combination of basis functions where N p is the number of basis under polynomial order N and c n is the modal basis while ' k i is the nodal basis. Piece-wise polynomial on each cell composite is the entire solution, namely The choice of basis function is critical to space discretization. Here, the discontinuous and continuous basis function space is defined as Following the method introduced in the work by Hesthaven and Warburton, 11 a set of hierarchical orthogonal basis is constructed by L 2 norm-based Gram-Schmidt method. For high-dimensional cases, expansion of basis function can be obtained by the tensor product of one-dimensional polynomial instead of generating specific function space.
Till now, we have an arbitrary order modal basis to represent the continuous or discontinuous function space. Although they are equivalent mathematically, the nodal basis is more favorable for space discretization for its convenience on handling boundary conditions. By providing the position of the nodal point, namely, Legendre-Gauss-Lobatto (LGL) quadrature points, the Lagrangian nodal basis can be derived from a modal basis. The transformation matrix between nodal and modal basis is called Vandermond matrix.
LGL points are essential to keep the matrix well conditioned. Since the basis functions talked above are defined on each cell piece-wisely, they maintain discontinuity on cell interfaces. However, to satisfy C 0 constraints in a continuous framework, the variable value of two adjacent cells should be the same on cell interfaces.

Variational form of discretization terms
Convective term. The convective term is in the momentum equation of the first substep. To increase stability for convection dominant flow, upwind DG method is used in the discretization of the convective term. Integrating by parts and using Lax-Friedrichs flux, the discrete formulation of the convective term can be derived as where l refers to maximum eigenvalue of respective flux Jacobian. In this case, l is approximated by the maximum velocity magnitude along flux direction on entire ∂D k . Boundary conditions are weekly imposed by given corresponding u + value on ∂O h , which is set using a mirror principle Laplacian term. The Laplacian term in pressure projection step is discretized with CG method. Since variables are continuous across interfaces in the framework of CG, flux integration can be ignored for internal faces.
The variational formulation is given as where the pressure gradient on Dirichlet boundaries is calculated by equation (8). Reference values on Neumann conditions are imposed by altering discretized linear system explicitly.
In the third viscosity substep, Laplacian term is discretized within DG framework using internal penalty formulation following methods used by Hesthaven and Warburton, 11 which helps to avoid spurious jumps on interfaces and provides sufficient stabilization. For the x component of velocity u h , the variational formulation of Laplacian term reads The penalty parameter t k for each element is defined as in which t is an empirical constant value and h is the mesh length scale of cell k. In the following context, t is set to a default value of 20 taken from work by Hesthaven and Warburton. 11 Gradient and divergence term. According to the work by Fehn et al., 9 the implementations of gradient and divergence term are critical to avoid small time instability. So, we shall illustrate the discretization of these two terms in detail, and show the differences between various methods.
The implementation of our DG-CG method is similar to the work by Botti and Di Pietro. 10 The gradient of pressure is calculated explicitly by derivative of basis function as shown in equation (21).
Divergence of intermediate velocity uses integration by parts to include information from neighbor cells and boundary conditions, as shown in equation (22). This is to make the divergence term consistent with viscosity substep (equation (19)) The approximation of velocity on cell interfaces are evaluated by central flux Imposing the boundary conditions illustrated in equation (17), the final formulation of divergence term is derived as where g u is the fixed value of velocity for Dirichlet boundaries.
To give a more direct illustration of the implementation, the corresponding pseudocode in Matlab style for two-dimwnsional (2D) case is provided as follows: In the pseudocode, Dr and Ds indicate the derivative matrixes in strong form, while Drw and Dsw are the derivative matrixes in week form. rx, sx, ry, and sy are the Jacobian coefficients. LIFT is the facial integration coefficients times inverse of mass matrix, Fscale is facial Jacobian coefficients, and variables with prefix fluxÀ are central flux for pressure or velocity. The calculation of pressure gradient uses only the information in local cell while velocity divergence needs one layer of neighbor cells.
To examine the small time stability of DG-CG method, it is necessary to find a reference. The method by Hesthaven and Warburton 11 is a widely used pure DG method, but reported to be unstable for small time step in the work by Krank et al. 5 In HW's method, both pressure gradient and velocity divergence are calculated explicitly, as shown in the following equation The corresponding pseudo-code is listed below.
As argued by Fehn et al., 9 the DG formulation of velocity divergence and pressure gradient term are the essential reasons for the small time instability. He recommends that integration by parts coupled with a suitable definition of boundary conditions should be used to keep stable and robust. The major difference 1 %compute gradient of pressure 2 pr = Dr*P; ps = Ds*P; 3 gradPx = rx.*pr + sx.*ps; 4 gradPy = ry.*pr + sy.*ps; 5 %compute divergence of velocity 6 ur = Drw*Ux; us = Dsw*Ux; 7 vr = Drw*Uy; vs = Dsw*Uy; 8 DivU = -(rx.*ur + sx.*us + ry.*vr + ... sy.*vs) + LIFT*(Fscale.*fluxU); 1 % compute gradient of pressure 2 pr = Dr*P; ps = Ds*P; 3 gradPx = rx.*pr + sx.*ps; 4 gradPy = ry.*pr + sy.*ps; 5 % compute divergence of velocity 6 ur = Dr*Ux; us = Ds*Ux; 7 vr = Dr*Uy; vs = Ds*Uy; 8 divU = rx.*ur + sx.*us + ry.*vr + sy.*vs; from HW's method is performing integration by parts for both pressure gradient and velocity divergence, which are consistent with schemes used in pressure projection and velocity correction substeps. In Fehn's method, the discretization of velocity divergence is almost the same as the one used in our DG-CG method (equation (22)). But the pressure gradient term is different, which is discretized with integration by parts. We can obtain the variational formulation as where g p is the fixed value of pressure gradient for Neumann boundaries. The pseudo-code is listed below.
It should be noted that the pressure in DG-CG method is continuous across element interfaces. That is to say, applying integration by parts for pressure gradient term in DG-CG method will not bring much difference. In the view of pressure gradient and velocity divergence, the DG-CG method is nearly the same as the one used in Fehn et al. 9 As Fehn claims that not applying integration by parts will lead to small time instability for pure DG method, one may ask whether the same phenomenon will happen in our mixed DG-CG method. Thus, we will also examine the DG-CG method using inconsistent velocity divergence scheme, as illustrated in equation (24). This method will be short for DG-CG_I in the later context.

Numerical results
The solver is implemented within the open-source framework HopeFOAM, 13 which aims at bringing high-order methods to official OpenFOAM based on second-order finite volume (Jasak et al., 14 Yang et al., 15 Wang et al. 16 ). The released HopeFOAM-0.1 provides high-order DGM, as well as high-order curved mesh, making it capable of achieving accurate space discretization. It inherits the user interfaces from OpenFOAM, including domain specified language, control dictionary, and boundary condition settings. Thus, developing high-order solver and cases within HopeFOAM is relatively convenient.
The linear system solver in HopeFOAM is based on PETSc, 17 which is a library intended for solving largescale scientific applications modeled by partial differential equations. It supports most of the commonly used linear solver and preconditioner for the sparse matrix. In the following context, the ilu preconditioned GMRES Krylov iterative method is used.

Pearson vortex
Pearson vortex simulation is an unsteady INS problem and frequently used to demonstrate the small time step instability. The solution is constituted by making temporal derivatives in momentum equation equal with viscous term and non-linear convective term balance with pressure gradient. The analytical solution is given by The computational domain is shown in Figure 1  to start from 2 À1 to 2 À16 decreased by a ratio of 0:5. To make a direct comparison, different methods discussed in the previous context are all performed on the test, including our DG-CG, DG-CG_I, the method used by Hesthaven and Warburton, 11 and the one proposed by Fehn et al. 9 The results in Figure 2 reproduce the phenomenon reported for Hesthaven and Fehn's method. For the HW's method, instabilities exist in all cases and the relative error increases dramatically at small time step size. For Fehn's method, integration by parts in pressure gradient and velocity divergence terms successfully remedy the instability for velocity. The relative error of velocity reduces to a stable value, which is dominated by spatial error.
The equal order DG-CG and DG-CG_I methods show great stability at small time step for both velocity and pressure at all polynomial orders. Similar relative L 2 error is obtained for these two methods for polynomial order 2 or 4. However, a significant increase of L 2 error of pressure is observed in Figure 2(b) for N = 1, while the increase in DG-CG is nearly not noticeable. This indicates that consistent velocity divergence scheme helps to obtain more accurate results for DG-CG method when the spatial resolution is extremely low. In our following context, the consistent velocity divergence scheme is used in our DG-CG method.
It should be noted that the L 2 error of pressure from Fehn's method is larger than DG-CG method in extremely low spatial resolution, namely, N = 1 and N = 2. Fehn pointed out that the increase of L 2 error is an indication of inf-sup instabilities. Although equal order DG-CG method has not been proved to satisfy the inf-sup condition, 18 its L 2 error is relatively smaller and bounded at small time step. Comparing to pure DG method, DG-CG method is observed to be much more accurate for the small time step.
Eigenvalue spectrum analysis. The eigenvalue spectrum analysis of the propagation matrix is performed to examine temporal stability of DG-CG method. Similar analysis has also been used in Leriche et al. 19 for spectral element collocation schemes and in Fehn et al. 9 for pure DG scheme. To perform this, the source term is discarded, harmonic boundary conditions are used for velocity, and the convective term is neglected. Then, the velocity U n + 1 can be derived by propagation matrix and U n , which reads To compute A, INS equations are treated as a whole system and the pressure is eliminated as an auxiliary variable represented by velocity. The stable criteria for this system is given as The max eigenvalues as a function of time step size are calculated on the problem defined in the previous subsection. Figure 3 shows the results of DG-CG method and pure DG by Hesthaven. For our DG-CG method, all eigenvalues are less than 1 and within the stable region. When the time step size comes to the limit Dt ! 0, the temporal derivative is dominating and max eigenvalue tends to 1. Thus, we can draw a conclusion that DG-CG method has formed a stable temporal  iterative system. On the contrary, the max eigenvalue of HW's pure DG method exceeds 1 at a small time step size, which indicates an unstable algebraic system. This is consistent with our findings from Figure 2.

Spatial stability in low viscosity situation
In this section, the spatial stability of equal order DG-CG method in low viscosity situation is analyzed. The behavior of equal-order DG-CG method at high Reynolds number is what we care about. In addition, how to keep stable is quite important when we try to apply the DG-CG method. This test is carried on Poiseuille flow case. In the Pearson vortex case, the viscosity term is set to balance with the convective term especially. Therefore, solving of pressure does not depend on viscosity. Thus, Pearson vortex is not suitable to investigate the instability caused by low viscosity. In Poiseuille flow, the viscous equals with pressure gradient term in the final steady state. Notably, the flow pattern of Poiseuille flow is common in computational fluid dynamics (CFD). This is a steady INS problem in two dimensions with an analytical solution The computational domain O = ½ÀL, L 2 with L = 1. The velocity profile is a hyperbolic equation with max value U max = 2m=s along the x axis. The inflow and up and down walls are set with Dirichlet boundary conditions, while outflow part is prescribed as Neumann boundary conditions. The flow starts from zero and ends at final time T = 30s, which ensures the flow has reached a steady state. The simulations are performed with equal order DG-CG method with consistent velocity divergence scheme, and the temporal scheme is set to J = 2.
First, the simulations are carried out with n = 5 3 10 À4 to examine the instability under high Reynolds number. Results are shown in Figures 4 and   1E

5, with various mesh scales and polynomial orders.
Numerical oscillations and unreasonable pressure fields can be seen in Figures 4(a) and 5(a), which indicate the instability of equal order DG-CG method. The instability is more significant under low spatial resolution configurations, namely, h = 1 and N = 2 in our results. According to the work by Lederer et al., 20 this instability caused by the velocity is not exactly divergence free.
In the situations of low viscosity, the error from velocity divergence is amplified leading to the unreasonable pressure field. Using nonconforming Crouzeix-Raviart element 21 or projecting the velocity to exactly divergence free space 6 could help to solve this instability. Besides, comparing Figures 4(c) and 5(c), we can find that finer mesh (h = 0:25) or higher polynomial order (N = 4) remedy these instabilities and greatly improve the numerical accuracy. Then, various time step sizes and dynamic viscosities are also taken into consideration. The relative L 2 error of pressure field at time t = 30s is collected in Table 1, covering h refinement, p refinement, various time step size, and viscosity. Simulations with n = 0:05 reach more accurate results comparing to cases with n = 5 3 10 À4 under the same mesh scale and polynomial order. This indicates that equal order DG-CG is more robust for relative low Reynolds number problem. And the high spatial resolution is required to avoid instability and obtain accurate results under this condition. Considering the effects of time step size, it should below the CFL limit to ensure convergence in the very beginning of the simulation, which is dominated by convection. Time step size dt = 10 À3 and dt = 10 À5 get very close L 2 error. However, slight increase can still be found in cases with dt = 10 À5 , which may relate to the inf-sup instability. When polynomial order N = 4, this increase seems more notable. Overall, the error introduced by a small time step size is bounded and do not affect the convergent. Otherwise, dt = 10 À5 is too small for the simulation, and setting the step size under CFL condition seems to be a good choice.
At last, we can conclude that equal order DG-CG method may fail for high Reynolds number cases under low spatial resolution. But with sufficient mesh resolution and polynomial order, coupling with proper time step size (CFL condition), accurate simulation results can still be achieved.

3D Beltrami flow and convergence test
This 3D case is developed by Ethier and Steinman 22 for benchmarking purposes, and is a counterpart of 2D Taylor vortex problem. This problem is defined on domain O = ½À1, 1 3 . The exact solution of velocity and pressure is given as u = À a e ax sin ay + dz ð Þ+ e az cos ax + dy ð Þ ½ e Ànd 2 t v = À a e ay sin az + dx ð Þ+ e ax cos ay + dz ð Þ ½ e Ànd 2 t w = À a e az sin ax + dy ð Þ+ e ay cos az + dx ð Þ ½ e Ànd 2 t p = À a 2 2 e 2ax + 2 sin ax + dy Here, the parameters are defined as dynamic viscosity n = 1:0, a = p=4, and d = p=2. To examine the spatial convergence, uniform structural meshes with 4 3 , 8 3 , and 16 3 elements are used. The point-wise error in L 2 norm of velocity and pressure in T = 1s are gathered in Figure 6. The slop between error and grid scale indicates the convergence order. Both velocity and pressure have reached theoretical convergence rate h N or even higher, compatible with pure DG discretization results mentioned in Hesthaven and Warburton. 11 To reduce the influence of error introduced by spatial discretization, the polynomial order is set with N = 5. A reverse Cuthill-McKee reorder method preconditioned LU-factor direct linear solver is chosen to solve the linear system. Temporal integration schemes with J = 1 and J = 2 are performed. Figure 7 shows the error in L 2 norm of velocity and pressure. For J = 1 case, the convergence rate is Dt, while J = 2 scheme reaches optimal convergence rate of Dt 2 . It should be noted that in the case of J = 2, the error does not continue to reduce when Dt\2E À4 since it is dominated by spacial accuracy.

Cavity flow
Here, we consider two-dimensional lid-driven cavity flow with relative high Reynolds number. The computational domain is a unit square. Upper bound (the lid) is enforced with a constant velocity u = (1, 0), while no-slip boundary conditions are set on other sides. The dynamic viscosity is configured with n = 10 À3 , yielding a Reynolds number for the flow of Re = 1000. Fifthorder DG-CG space discretization and J 2 time integration scheme are applied on a triangular mesh with 3606 cells. The spatial resolution is sufficient to avoid low viscosity related instability and obtain accurate results. The simulation advances until reaching a steady state. Figure 8(a) shows the streamline of steady-state cavity. Two vortices appear at the bottom corner, which is consistent with the literature. The profile of x and y components of velocity along vertical and horizontal center line is drawn in Figure 8(b). The simulation results (the lines) fit well with the reference value (dots). 23 With special care on the mesh resolution and discretization order, DG-CG method can simulate relative high Reynolds problem and keep stable.

3D laminar flow past cylinder
This is a very challenging case with unsteady timevarying inflow profile. As a simplification of practical applications, lots of papers have used this case to verify the 3D simulation accuracy of CFD code. For the benchmark problem, a cylinder locates at the middle line of a channel, while no-slip boundary conditions are employed on the walls. The detailed geometry configuration can refer to Scha¨fer et al. 24 and John 25 The inlet velocity is set with where U m = 2:25m=s. The simulation starts with zero inflow at t = 0 and finishes at t = 8s. Guided by the findings in ''Pearson vortex'' and ''Spatial stability in low viscosity situation'' sections, the simulation is successfully carried out by DG-CG method with up to fifth-polynomial order. Although using a coarse mesh composed of 50,000 tetrahedral, results are nearly converged at third order due to the high-order method. Special attentions are paid on the cylinder boundary, where high-order curved elements are employed on the nearest cells. Figure 9 depicts the velocity field at t = 2:4s using DG-CG method with polynomial order N = 3, which is indistinguishable from results with N = 4, 5.
The maximum drag coefficient and minimum lift coefficient are major comparing criteria for this case. Results obtained by DG-CG and pure DG method under different polynomial orders are listed in Table 2. Accurate results are gained under both methods according to the reference values in Bayraktar et al. 26 It is worth to mention that results by our high order methods with N = 3 are as accurate as OpenFOAM (second-order FVM) on a mesh with 3 million cells. DG-CG method uses about 33% degree of freedom (Dof) than FVM to reach similar accuracy. Performance evaluation DG method is preferred for its high-order convergence rate, stability, and high parallel efficiency while coming along with the disadvantages of extraordinary computation and memory overhead. This has been a major factor in obstructing its spread and application on complex problems. In the splitting method, solving of pressure takes up the major part of computational overhead. The solver with discontinuous velocity and continuous pressure has been examined to achieve high-order spacial convergence rate in previous sections. Moreover, pressure discrete by CG has the potential to be more efficient comparing to the discontinuous method. In this section, the performance of equal order DG-CG and pure DG method are analyzed under ilu preconditioned GMRES solver with the Pearson vortex case. The test ends at final time T = 1:0s with time step set by CFL condition.
Although it is a very simple case, results can indicate the general tendency and provide references for other cases.
The performance statistics of DG-CG and pure DG method are collected in Table 3, including the iterative number (n_it), computational time, and the number of degree of freedom (nDof) for each substep. Notation ''pressure-cg'' indicates the substep of solving pressure in DG-CG method, ''pressure-dg'' is a similar substep in pure DG method, and ''velocity-dg'' is the substep of solving velocity correction in both methods. Similar to the LU direct solver cases, DG-CG method has fewer Dofs and need less iterative steps to get convergence comparing to pure DG method. The speedup reaches up to 4 times at 896 test case while 6 times at 3584 test case, which is more significant at high-order scenario.  Figure 9. Snapshot of velocity magnitude for 3D laminar flow past cylinder using DG-CG method with J = 2 and polynomial order N = 3 at t = 2:4s.

Conclusion
In this article, we analyzed in detail the stability of the splitting method with discontinuous velocity and continuous pressure. Results are also compared to pure DG method by Hesthaven and Warburton 11 and Fehn et al. 9 For the small time step size problem, we first reproduce the instability with HW's method on the simulation of Pearson vortex, while our DG-CG method remains stable at the same settings. As the time step size decreases continuously, the relative L 2 error of pressure increase litter. In addition, the consistent velocity divergence discretization scheme helps to get more accurate results, especially for low polynomial order cases. Second, we demonstrate the instabilities of equal order DG-CG method in high Reynolds number cases. Through comprehensive experiments of Poiseuille flow covering various viscosities, time step sizes, mesh scales, and polynomial orders, we numerically show that sufficient polynomial order and mesh helps remedy the instability and obtain accurate results for high Reynolds flow.
In the work of Botti and Di Pietro, 10 the accuracy and performance of high-order polynomial are not discussed, where the highest pair is DG(3)-CG(2). In the stability analysis, we numerically show that higher polynomial order can lead to a more robust scheme in low viscosity situations. By the example of 3D Beltrami flow, up to sixth-order spatial convergence rate and second-order temporal convergence rate have been demonstrated. The cavity case is run on relative high Reynolds number flow condition, which confirms the stability of DG-CG method. At last, we check out the performance of DG-CG method at high polynomial order comparing with pure DG. Results show that DG-CG method is more efficient than pure DG method, to be specific, a speedup of 2.5 to 6 times on ilu preconditioned GMRES solver, which meets our initial expectation to develop mixed DG-CG solver.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.