Adaptive local discontinuous Galerkin methods with semi-implicit time discretizations for the Navier-Stokes equations

In this paper, we present a mesh adaptation algorithm for the unsteady compressible Navier-Stokes equations under the framework of local discontinuous Galerkin methods coupled with implicit-explicit Runge-Kutta or spectral deferred correction time discretization methods. In both of the two high order semi-implicit time integration methods, the convective flux is treated explicitly and the viscous and heat fluxes are treated implicitly. The remarkable benefits of such semi-implicit temporal discretizations are that they can not only overcome the stringent time step restriction compared with time explicit methods, but also avoid the construction of the large Jacobian matrix as is done for fully implicit methods, thus are relatively easy to implement. To save computing time as well as capture the flow structures of interest accurately, a local mesh refinement (h-adaptive) technique, in which we present detailed criteria for selecting candidate elements and complete strategies to refine and coarsen them, is also applied for the Navier-Stokes equations. Numerical experiments are provided to illustrate the high order accuracy, efficiency and capabilities of the semi-implicit schemes in combination with adaptive local discontinuous Galerkin methods for the Navier-Stokes equations.

was firstly introduced in 1973 by Reed and Hill [1] in the framework of linear neutron transport and then further extended by Cockburn et al. [2][3][4][5] for solving nonlinear hyperbolic conservation laws through a combination of explicit, Strong-Stability-Preserving Runge-Kutta (SSP-RK) time discretizations and total variation bounded nonlinear limiters to achieve non-oscillatory properties for strong shocks. The DG methods combine two advantageous features commonly associated to continuous finite element and finite volume methods, such as high order accuracy, flexibility on complex geometry, flexibility for h-p adaptivity, and so on, and are indeed a natural consideration when solving hyperbolic conservation laws, such as the compressible Euler equations. However, when it comes to partial differential equations (PDEs) containing high order spatial derivatives, e.g., the compressible NS equations where the viscous and heat fluxes exist, a severe difficulty with the approximations of the numerical fluxes for solution derivatives arises by the direct application of the DG methods, and a naive arithmetic mean of the solution derivatives without taking into account the possible jump of the solutions will yield a weakly unstable scheme.
In order to properly resolve the solution derivatives in the NS equations at the interfaces, plenty of numerical methods have been proposed in the literature. In 1997, Bassi and Rebay [6] firstly attempted to apply the DG methods to the compressible NS equations, later on, they introduced an improved method [7] based on the previous scheme in order to maintain the compactness and stability for the pure diffusion problems. Hartmann and Houston [8] proposed employing the symmetric interior penalty method, which can guarantee optimal order of convergence in terms of L 2 norm of the error, for the discretization of the leading order terms of the compressible NS equations. To deal with moving and deforming boundaries, Klaij and van der Vegt [9] presented a space-time DG method for the compressible NS equations. The key feature of the spacetime DG method is that no distinction is made between space and time variables and thus this provides optimal flexibility to deal with time dependent boundaries and deforming elements. Based on a so-called inter-cell reconstruction, Luo et al. [10] developed a reconstruction-based DG method for the NS equations on arbitrary grids. The numerical viscous and heat fluxes at interfaces are obtained by locally reconstructing a smooth solution with a least-square method from the underlying discontinuous solution. Based on a direct weak formulation for the parabolic equations, a so-called direct DG method for solving second-order diffusion problems was originally introduced by Liu and Yan [11], and further extended by Cheng et al. [12] to solve the compressible NS equations.
Inspired and attracted by the high order accuracy, easy extension to high order PDEs and arbitrary grids of the LDG method, we select it as the spatial discretization in our study. Besides, adopting h-adaptive technique and semi-implicit time marching methods are another two aspects we focus on in this paper to improve the computational accuracy and efficiency as well as deduce the computational complexity. The LDG method was initially developed to solve nonlinear convection-diffusion equations by Cockburn and Shu [13], enlightened by Bassi and Rebay's successful simulation for the compressible NS equations. The idea behind the LDG method is to properly rewrite the PDEs containing high order derivatives as a first-order system by introducing auxiliary variables, then apply the DG method to the system, in which correctly designing the numerical fluxes is the key ingredient to guarantee stability and local solvability of the auxiliary variables. There has been abundant literature on designing and analyzing the LDG schemes for different types of high order PDEs, for a detailed review, see [14] and the references therein.
By sharing the advantages of the DG method, the LDG scheme requires no imposition of solution continuity between adjacent elements and allows the appearance of hanging nodes, thus is extremely flexible to do local grid refinement. To start the mesh adaptation, a criterion to determine whether an element in the computational mesh needs to be refined or coarsened is demanded. Generally, the criteria can be divided into two types, i.e., error estimators and heuristic indicators [15]. Error estimators are based on a so-called a posteriori error estimation which seeks to bound the error with respect to a given norm. As a comparison, heuristic indicators usually depend on the local gradient of thermodynamical variables such as the density, pressure, entropy and so on, or the local divergence or curl of the velocity field, thus are relatively simple to be implemented for complicated PDEs. In the current, a rigorous mathematical analysis to develop the a posteriori error estimation for the LDG scheme of the compressible NS equations is out of our scope, therefore we consult the criteria in [16,17] for the compressible Euler equations and choose a combination of different heuristic indicators. In practical implementation, the candidate elements marked for refinement are indeed to be refined, however, to retain the mesh quality, some other elements around these elements may also require to be refined. In contrast, to achieve the same goal, the candidate elements marked for coarsening may not be coarsened depending on their neighbors. What's more, to maintain conservation and accuracy of the numerical solutions during refinement and coarsening, the L 2 projection is implemented for the prolongation and restriction of grid.
The application of any above extension of the DG methods to the compressible NS equations will generate a large coupled system of ordinary differential equations (ODEs), which require accurate and efficient time integration methods to march in time. The explicit, high order SSP-RK time discretizations, which are suitable choices for hyperbolic conservation laws, will sustain extremely small time step restriction for stability, but not for accuracy due to the appearance of stiffness terms of the NS equations, i.e., the viscous and heat fluxes terms. Nearly all the papers referred above applied the implicit temporal discretizations, e.g., the backward Euler method, to overcome this restriction. This kind of remedy is efficient to some extent, since the Courant-Friedrichs-Lewy (CFL) number can be adjusted to be quite large, but the construction of the large Jacobian matrix and the choice of efficient preconditioners for the linear systems arising in the inner loop greatly increase the difficulty in implementation. In addition, for unsteady flow problems with complicated flow structures, a time step of at least the same order of magnitude as the mesh size is required to capture the complex flow structures accurately at any time. By turning to the semi-implicit time discretization methods, specifically, IMEX-RK and SDC methods, and treating the convective, viscous and heat fluxes differently, no Jacobian matrix is required to be assembled and only two linear systems for the momentum and energy equations, respectively, need to be solved. Numerical experiments show that the residual of the GMRES solver even without any preconditioner can approach the machine error within a dozen of, sometimes two or three, outer iterations if the mesh quality is not too bad. L. Pareschi and G. Russo considered IMEX-RK methods in time for hyperbolic systems of conservation laws with stiff relaxation terms in [18]. In [19,20] Wang et al. analyzed the stability and error estimates of the LDG methods coupled with IMEX-RK time discretizations for the linear convection-diffusion equations in one and two dimen-sions, and unconditional stability and optimal error estimates are obtained. Xia et al. [21] explored three different time discretization techniques including SDC methods for solving the stiff ODEs resulting from an LDG spatial discretization to PDEs containing high order spatial derivatives, and verified that all the three methods are efficient and a time step t = O( x) is allowed for k-th order PDEs. The LDG method is more flexible for treating the nonlinear terms when introducing the auxiliary for high order derivative terms. Due to the local properties of the LDG methods, the resulting implicit scheme is easy to implement and can be solved in an explicit way when it is coupled with iterative methods, which has been demonstrated in the references [22,23]. For other work on semi-implicit schemes for gasdynamics, we suggest the readers consulting [24][25][26][27][28][29][30] and the references therein.
The outline of this paper is organized as follows. Governing equations are given in Section 2. In Section 3, we present the LDG scheme for the compressible NS equations, including the detailed and clear treatment for the numerical fluxes for subsonic inflow/outflow, supersonic inflow/outflow and solid wall boundary conditions. Section 4 is devoted to the discussion about two kinds of semi-implicit time discretization methods we adopt in this paper, i.e., IMEX-RK and SDC methods, and the fully discrete schemes of the NS equations. In Section 5 we provide the implementation details of the h-adaptive technique. Section 6 contains numerical results for the unsteady compressible NS equations, which demonstrate the high order accuracy and capabilities of the presented methods. Finally we give conclusions in Section 7.

Governing equations
Let ⊂ R d be a bounded domain with d ≤ 3, the NS equations governing the dynamics of viscous compressible flows express the conservation of mass, momentum and energy, and can be written in dimensionless and conservative form [31] ∂U ∂t in ×(0, T], T > 0, with U ∈ R d+2 the vector of conserved variables and F, G ∈ R (d+2)×d the inviscid and viscous fluxes, which are defined as respectively. Here ρ, u ∈ R d and p denote the mass density, velocity field and pressure, respectively. I represents unit tensor. The total energy per unit volume E is the sum of internal and kinetic energy with e the specific internal energy. In order to close the system, we consider the equation of state for a calorically perfect gas c v is the ratio of c p and c v , the specific heats at constant pressure and volume, respectively. By Newtonian approximation and Stokes hypothesis, the viscous stress tensor τ ∈ R d×d relating to the derivatives of the velocity field u is defined as where the dynamic viscosity coefficient μ, which is determined through Sutherland's law, reads in dimensionless form with Re ∞ the Reynolds number, T the temperature, T s a constant and subscript ∞ denoting the uniform free-stream values. The heat flux vector q ∈ R d caused by the gradient of temperature T is defined as Pr is the thermal conductivity coefficient, with Pr the Prandtl number. In addition, the specific internal energy e can relate to the temperature T by e = c v T.
As a matter of fact, the specific heats c p and c v in dimensionless form can be defined as with Mach number M ∞ for uniform flow given by In this paper, unless otherwise stated, we will use γ = 1.4, Pr = 0.72 for air, and μ as a constant for simplicity, i.e., μ = 1 Re ∞ .

The LDG method
In this section, we will introduce the LDG method for the spatial discretization of the compressible NS equations. Special attention will be paid to the numerical fluxes on various boundaries.

Notations and finite element spaces
For the DG spatial discretizations, the domain is approximated with a tessellation T h consisting of nonoverlapping shape-regular elements K, which satisfy the condition K∈T h K := h → as h → 0, with the mesh size h := max K∈T h h K and the diameter h K of element K. Let denote the union of the boundaries of all the elements K and classify them into internal faces I and boundary faces B , i.e., = I ∪ B . We also assume that B may be decomposed as follows where sub−in , sub−out , sup−in , sup−out and W are distinct subsets of B representing subsonic-inflow, subsonic-outflow, supersonic-inflow, supersonic-outflow and solid wall boundaries, respectively. For solid wall boundaries, we also distinguish them either according to slip (reflective) and no-slip conditions, i.e., W = W ,slip ∪ W ,no−slip , or according to isothermal and adiabatic conditions, i.e., Let e ∈ I be an internal face shared by the "left" and "right" elements K L and K R , i.e., e = ∂K L ∩ ∂K R , where the so-called "left" and "right" can be uniquely defined for each internal face according to any fixed rule. Suppose φ is a function on K L and K R , but possibly discontinuous across e, let φ L and φ R denote (φ| K L )| e and (φ| K R )| e , the left and right traces, respectively. Similar definitions can be obtained component-wisely for vector-valued and matrix-valued functions.
For the LDG discretizations, we require to introduce the finite element spaces. Actually each element K of the tessellation T h is connected to a reference elementK through some mapping F K . The mapping F K :K → K from the reference elementK to the real physical element K is a function defined in the space of the reference element for each independent variable. For example, for a quadrilateral element K in two dimension,K =[ −1, 1] 2 is the unit square and F K is expressed in terms of the nodal shape functions N i , i = 1, · · · , 4 by are the four vertex coordinates of K and see schematic map in Fig. 1. Then the finite element spaces associated with the tessellation T h are given by with P k (K) the space of polynomials of degree up to k with respect to d variables in the reference elementK. Note that the functions in V h , W h and h are allowed to be completely discontinuous across element interfaces. Further, we also define the inner product notations in element and on face as for the scalar-valued functions φ, ϕ, vector-valued functions u, v and matrix-valued functions σ , η, respectively, where

The LDG discretization
In order to propose the LDG discretization for the NS equations, we firstly rewrite (1) as a first order system, which is composed of the primary equations ∂U ∂t and the auxiliary equations with tr(z) := z ii the trace of the matrix z. Based on the weak formulation of (2)− (3), which is obtained by multiplying (2)−(3) by test functions, integrating over some domain, and then performing an integration by parts, we can obtain the following semi-discrete LDG scheme : find U h ∈ V h , z h ∈ h and q h ∈ W h , such that for all test functions V ∈ V h , σ ∈ h and p ∈ W h , the following equations are satisfied and n is the unit outward normal vector to the boundary ∂K. The terms denoted with a hat in (5) and (6) in the cell boundary terms from integration by parts are the so-called numerical fluxes, which are functions defined on the faces and should be designed based on different guiding principles for different PDEs to ensure stability and local solvability of the intermediate variables. Consider a face e ∈ ∂K, and we denote by superscripts L and R the internal interface state and neighboring element interface state, respectively. Note, the neighboring element could be a ghost element which lies exterior to the computational domain.
• If e ∈ I is an internal face, we choose the local Lax-Friedrichs flux for the convective part and central flux for the other flux terms. In detail, where α e is the biggest eigenvalue of the Jacobian matrix ∂(F h n) ∂U on e. • If e ∈ sub−in ∪ sub−out ∪ sup−in ∪ sup−out is a farfield boundary face, we define the right external states as then the numerical flux F is similarly defined as in the internal face, and is a solid wall boundary face, in order to weakly prescribe the slip and no-slip boundary conditions, i.e., u · n = 0 and u = 0, respectively, for the convective flux F, we again define the right external states as then F can be formulated the same as in the internal face. The velocity at the boundary is defined as In addition, for the isothermal case e ∈ W ,iso , T and G are given by with T w the given temperature on the wall. For the adiabatic case e ∈ W ,adia , T and G are given by Now we have completed the definition of the LDG scheme for the NS equations.

Time discretizations
The LDG spatial discretization for the NS equations typically results in a system of ODEs which contains a non-stiff term F N (t, u) as well as a stiff term F S (t, u). One would then need to use a suitable ODE solver to discretize the temporal variable. Generally, in many cases, the non-stiff component F N (t, u) is nonlinear and the stiff component F S (t, u) is linear, such as Burgers equation with viscous term, KdV type equations, etc. It would be desirable to treat the different components separately, more specifically, to treat the nonstiff term explicitly and stiff term implicitly. This kind of treatment can not only relax the severe time step restriction due to implicit integration for the stiff term, but also be easy to implement since we apply explicit discretization for potentially nonlinear nonstiff terms. And what we would like to emphasize is that, even though the stiff term of the NS equations is nonlinear with respect to all the conserved variables, we can still update the variables at the new time level by solving a series of linear systems due to the special structure of the NS equations. In the following work, we will consider two kinds of time discretizations: IMEX-RK methods and semi-implicit SDC methods.
To numerically solve (7), we divide the time interval (0, T] into M elements by the partition 0 = t 0 < t 1 < · · · < t n < · · · < t M = T and denote by t n = t n+1 − t n the time step at level n and u n the numerical approximation of u(t n ).

IMEX-RK discretizations
By applying the IMEX-RK time marching methods, the numerical solution of (7) advanced from time t n to t n+1 is given by can interpret the IMEX-RK methods simply and clearly using a Butcher tableau For a detailed introduction to IMEX-RK schemes, we refer the readers to [18,32,33]. The followings are examples of the first, second and third order IMEX-RK methods, respectively: • First order (one stage): It is obvious that the first order IMEX-RK method is just taking the forward Euler discretization for the non-stiff term and the backward Euler discretization for the stiff term.

Semi-discrete SDC discretizations
Algorithm 1 Semi-discrete SDC method Compute the initial first order approximation u 1 n,0 = u n for m = 0, · · · , P − 1 do u 1 n,m+1 = u 1 n,m + t n,m F N t n,m , u 1 n,m + F S t n,m+1 , u 1 n,m+1 end for Compute successive corrections The SDC time discretization is based on the Picard integral equation and low order time integration methods, which are corrected iteratively, with the order of accuracy increased by one for each additional iteration. The key advantage of SDC method is that it can systematically and simply develop time integration methods of arbitrary high order of accuracy.
To define the semi-discrete SDC scheme, we further divide the time interval [ t n , t n+1 ] into P subintervals by points t n,m for m = 0, 1, · · · , P such that t n = t n,0 < t n,1 < · · · < t n,m < · · · < t n,P = t n+1 . Let t n,m = t n,m+1 − t n,m and u k n,m denote the k-th order approximation to u(t n,m ). We choose the points {t n,m } p m=0 as Gauss-Lobatto nodes on [ t n , t n+1 ]. Then the algorithm to calculate (k + 1)-th order accuracy numerical solution of (7) advanced from time t n to t n+1 is given in Algorithm 1, where I m+1 m F N t, u k + F S t, u k is the integral of the P-th degree interpolating polynomial on the P + 1 points t n,m , F N t n,m , u k n,m + F S t n,m , u k n,m P m=0 over the subinterval [ t n,m , t n,m+1 ].

Fully discrete schemes
Finally, as the implementation of IMEX-RK and SDC time discretizations to the semidiscrete LDG scheme (5) and (6), we will present the fully discrete LDG schemes in the following for the NS equations. We only present the third-order fully discrete scheme here and the lower order fully discrete schemes can be deduced similarly.
• Third order IMEX-RK-LDG scheme The LDG scheme with the third order IMEX-RK time marching method is given as : and notations * n,0 = * n .

Mesh adaptation
In this section, we mainly present the detailed implementation of the h-adaptation algorithm for the structured mesh consisting of quadrilateral elements in 2D or hexahedral elements in 3D, including the method to refine and coarsen elements, the criteria to select candidate elements for refinement and coarsening, and the specific procedure of mesh adaptation.

Refinement and coarsening of elements
For the data structure of adaptive mesh, we adopt a generalized binary tree, i.e., quadtree in 2D and octree in 3D, to store all the elements, that is, a parent element is always divided into 4 in 2D or 8 in 3D child elements. The DG method is extremely local in data communication and allows the appearance of hanging nodes. This means that an element can be refined an unlimited number of times, regardless of its immediate neighbors. However, if a huge difference between the levels of adjacent elements exists, the stability of the resulting numerical scheme will be damaged. In order to enhance the mesh quality and the stability of the scheme, we impose the quadtree or octree balanced, that is, the levels between two adjacent elements differ at most by 1.
The numerical solution U h restricted in element K can be expressed as where U K l (t) and v K l (x), l = 1, · · · , N b denote the degrees of freedom and basis functions, respectively, in element K. To maintain the accuracy and local conservation of the solutions, we adopt L 2 projection to obtain the degrees of freedom in the new generated elements during refinement and coarsening. In detail, provided the numerical solution U h is already known on the mesh T h (t n ), we require to determine the degrees of freedom U K l (t n ), l = 1, · · · , N b in the new element K ∈ T h (t n+1 ). Let U h be the L 2 projection of U h , which is computed through the following equations: Then the degrees of freedom U K l (t n ), l = 1, · · · , N b can be obtained by substituting (14) into (15).

Indicators
A criterion will be presented to initially determine the candidate elements for refinement and coarsening in a given mesh. According to [16,17], the gradient of density finds shock and contact discontinuity well, the divergence of velocity is direction independent and very effective in locating shock including strong shock and weak shock and the curl of velocity is also direction independent and very effective in finding shear and vortex. For these three different indicators, we compute the following quantities respectively: η gi = |∇ρ|d 3 2 i , η di = |∇ · u|d 3 2 i , η ci = |∇ × u|d 3 2 i , i = 1, · · · , N c , where N c is the total number of elements, d i = √ |K| or d i = 3 √ |K| depending on the dimension. The standard deviation of the gradient of density, the divergence and the curl of velocity are respectively. In our work, a single indicator or a combination of two of the above indicators will be taken into account for different problems. Suppose the level of all the elements in the initial mesh equals zero, and each element in the initial mesh can be refined at most LEV times. Based on the values of indicators, if only a single indicator, e.g., η gi is considered, then • If η gi > ω 1 η g and the level of K < LEV , then K is marked as a candidate element for refinement, • If η gi < ω 2 η g and the level of K > 0, then K is marked as a candidate element for coarsening, with problem dependent parameters ω l , l = 1, 2. If a combination of two indicators, e.g., η gi and η di are employed, then • If η gi > ω 1 η g or η di > ω 3 η d , and the level of K < LEV , then K is marked as a candidate element for refinement, • If η gi < ω 2 η g and η di < ω 4 η d , and the level of K > 0, then K is marked as a candidate element for coarsening, with problem dependent parameters ω l , l = 1, 2, 3, 4. In general, ω 1 and ω 3 are chosen between 1.1 and 1.5, ω 2 and ω 4 are chosen between 0.1 and 0.5.

Strategy for refinement and coarsening
In the above subsection, several candidate elements for refinement and coarsening are selected through given criteria. And whether these elements can be ultimately refined or coarsened, they also have to adhere to the principle that the level difference between two adjacent elements is at most 1 to improve the stability and accuracy of the LDG methods. To satisfy this principle, we adopt the strategy "refinement must, coarsening can" in [15], which means that a candidate element flagged for refinement is certainly to be refined, and in comparison, a candidate element marked for coarsening may not be coarsened.
After the execution of Algorithm 2 and Algorithm 3 in sequence, we can obtain two sets C r and C c which contain the final elements for refinement and coarsening, respectively. In Algorithm 2, all the candidate elements for refinement will be refined, and some neighbors of these elements are also added to the set C r to ensure the mesh quality. In Algorithm 3, a candidate element for coarsening may be deleted from the set C c due to either the absence of marked for coarsening for its brother elements, or the level difference with its neighbors.

Flow chart for mesh adaptation
As an end of this section, we present the general procedure for the adaptive LDG methods. Given the initial coarse mesh T h (t 0 ) and final computational time T,

Algorithm 2 Refinement of candidate elements
Suppose the set of all the candidate elements for refinement is C r for each element K in C r do for each neighborK of K do ifK / ∈ C r and the level ofK < the level of K then addK to C r end if end for end for

Algorithm 3 Coarsening of candidate elements
Suppose the set of all the candidate elements for coarsening is C c for each element K in C c do Set K p as the parent element of K if all the child elements of K p belong to C c then if the level of all the neighbors of K p ≤ the level of K then Step 1 Initialize the level of all the elements in T h (t 0 ) to be 0 and obtain the degrees of freedom in each element U K l (t 0 ), l = 1, · · · , N b , K ∈ T h (t 0 ) through L 2 projection for the initial exact solution.
Step 2 Given the mesh T h (t n ) and the degrees of freedom in each element is gained as follows: -Identify the candidate elements for refinement and coarsening in T h (t n ) through criteria in subsection 5.2. -Determine the final element sets C r and C c for refinement and coarsening through Algorithm 2 and Algorithm 3, respectively and subsequently. -For each element in C r , divide it into 4 in 2D or 8 in 3D child elements, get the degrees of freedom of child elements by use of L 2 projection, and increase the level of child elements by one. -For child elements with the same parent element in C c , delete them and get the degrees of freedom of the parent element through L 2 projection.
Step 3 Obtain U K l (t n + 1), l = 1, · · · , N b , K ∈ T h (t n + 1) by aid of the LDG method coupled with an IMEX-RK or SDC method.

Numerical experiments
In this section, we will carry out several numerical experiments including accuracy tests and other benchmark problems for the NS equations to evaluate the performance of the h-adaptive LDG methods. In addition, IMEX-RK and/or SDC time integration methods are adopted to solve the ODEs resulting from the spatial discretizations. Note, at every stage for the evolution in time, GMRES solver without any preconditioner is adopted to solve the two successive linear systems with respect to momentum and energy equations. In all experiments, the time step t n at time t n is chosen to satisfy the following CFL condition d t n min K h K max where d ≤ 3 denotes the dimension, h K is the diameter of element K, u K and a K are the velocity vector and speed of sound at element K, respectively, cfl is the CFL number taken as 0.98, 0.3, 0.18 for the first, second and third order time discretizations, respectively. In the following simulations, unless otherwise stated, an element in a mesh will be refined at most LEV = 3 times, that is, the depth of the generalized binary tree is 3.

Example 1 Accuracy test in two dimension
In order to demonstrate the uniform high order accuracy both in space and time of our proposed h-adaptive LDG scheme in conjunction with IMEX-RK or SDC methods, for the first test case we consider a smooth exact solution ρ(x, t) = 0.6 + 0.1 sin(5πt) cos(2πx) cos(2πy), u(x, t) = sin(3πt) sin(2πx) sin(2πy), v(x, t) = sin(3πt) sin(4πx) sin(4πy), (16) p(x, t) = 0.8 + 0.1 sin(πt) sin(2πx) cos(2πy), for the compressible NS Eqs. (1) with d = 2 and an additional source term, which is acquired by inserting (16) into (1). Additionally, we take the computational domain as [ 0, 1] 2 with periodic boundary conditions and Reynolds number Re ∞ = 200, 1000, 5000, respectively. All the simulations are performed until time T = 0.1 on a set of successively globally refined background meshes containing N × N, N = 16, 32, 64, 128 quadrilateral elements, respectively. For a fixed N × N background mesh, we randomly choose 2N candidate elements marked for refinement and the remaining elements marked for coarsening at each time step, and gain the final element sets C r and C c through Algorithms 2 and 3. The initial 16 × 16 background mesh and the adaptive mesh at some time based on the 32 × 32 background mesh are shown in Fig. 2. We adopt (k + 1)-th order IMEX-RK or SDC time marching methods when k-th degree polynomials are employed. The L 2 errors and orders of accuracy of density, velocity and pressure for both two schemes are presented in Tables 1 and 2, respectively, which illustrate that our two schemes can both obtain the expected optimal order of accuracy for these physical quantities with various Reynolds numbers. Moreover, by comparing the results in Tables 1 and 2, we can conclude that the IMEX-RK and SDC methods with the same order share nearly identical accuracy based on the background meshes with the same number of elements. For further efficiency comparison between IMEX-RK and SDC methods, we run the simulations on the uniform 64 × 64 mesh to eliminate the influence of frequent mesh refinement and coarsening. It can be seen from Table 3 that, compared with the IMEX-RK scheme, the SDC time discretization takes nearly the same time for the second order case and twice the time for the third order case, respectively. This can evidently be explained by their respective number of stages (there are 2 stages for both the second order IMEX-RK and SDC schemes, 3 stages for the third order IMEX-RK scheme and 6 stages for the third  order SDC scheme). However, unlike the IMEX-RK scheme, the SDC method of arbitrary high order of accuracy can be systematically and simply developed.

Example 2 Accuracy test in three dimension
In this example, we extend our accuracy test to three dimension, which is much more challenging to implement the procedure of refinement and coarsening for a general hexahedral mesh. To this end, similar to the two dimensional counterpart, an artificial exact smooth solution w(x, t) = e −2t sin(2πx) sin(2πy) cos(2πz), p(x, t) = 1.5 + 0.2e −t cos(2πx) cos(2πy) cos(2πz),  Fig. 3 for the 3D adaptive sample mesh. We also compute the L 2 errors and orders of density, velocity and pressure for IMEX-RK-LDG and SDC-LDG schemes. As in the previous case, Tables 4 and 5 show that our adaptive methods can also achieve optimal order of accuracy for 3D NS equations and nearly the same accuracy is obtained for the two schemes with identical setup.

Example 3 Efficiency test of h-adaptive methods
For the purpose of illustrating the superiority of the presented mesh adaptation technique compared with the uniform mesh method in terms of storage and computational time, we consider the evolution of an isentropic vortex in inviscid flows. This procedure is governed by the 2D Euler equations, which are obtained by neglecting the viscous effects and heat conduction, that is, omitting the right-hand side of Eqs. (1). The initial solution of this problem is obtained by adding some perturbations to the uniform mean flow The vortex perturbations of the temperature T, velocity field (u, v) and entropy S can be formulated as follows: is the coordinate of the vortex center and is the vortex strength. Then the initial solution is determined through isentropic relation The exact solution to this problem is simply the passive convection of the initial solution. We take (x 0 , y 0 ) = (10, 10), = 5 and computational domain as [ 0, 50] 2 . The numerical solutions are computed using piecewise quadratic polynomials and third order SSP-RK method which is more suitable for hyperbolic conservation laws. For more details about SSP-RK methods, see [2][3][4][5].  Table 6 quantitatively shows that our h-adaptive methods can save 89.9% in storage and 70.4% in CPU time compared with the uniform mesh method when they achieve the same L 2 error 4.84E-6 of density.  To explicitly show the shock-capturing property of our proposed methods, we simulate the Riemann problems which are usually taken as Euler benchmark problems directly using the NS equations. Specifically, the NS Eqs. (1) with d = 1 and initial conditions

Example 4 Shock-tube problems
are considered in this test case. We adopt the classical Sod and Lax initial settings: We run the two simulations with Reynolds number Re ∞ = 1800, P 2 element, third order IMEX-RK methods as well as a combination of indicators η gi and η di with parameters ω 1 = 1.2, ω 2 = 0.3, ω 3 = 1.2, and ω 4 = 0.3 on a background mesh of 200 cells till  time T = 0.25 and T = 0.13, respectively. Figures 6 and 7 show the comparisons of density, velocity, pressure, specific internal energy between the Euler exact solutions and the NS numerical solutions. Note, the points for the NS numerical solutions are plotted every four cells and more points are clustered near the shock waves, contact discontinuities and rarefaction waves due to mesh adaptation to get better resolution.

Example 5 Taylor-Green vortex
We consider the Taylor-Green vortex, which is a low Mach number problem, in this test case. The initial solution of this problem is with p 0 = 10 5 . This configuration leads to a characteristic Mach number M ≈ 1.7 · 10 −3 , thus corresponds to the low Mach number regime [27]. The computational domain is [ 0, 2π] 2 . We simulate this procedure on a 64×64 background mesh with periodic boundary conditions till final time T = 0.1. The (k + 1)-th order IMEX-RK time discretization will be employed when P k element is used. A single adaptive indicator η ci with parameters ω 1 = 1.5, ω 2 = 0.5 is adopted. As is shown in Fig. 8, the vortex can be accurately captured by the adaptive mesh, and the contours become smoother when higher order elements are employed.

Example 6 Subsonic flow around a NACA0012 airfoil
In this test case we consider a subsonic flow around a NACA0012 airfoil which is symmetric about the x-axis with upper and lower surfaces described by functions [34] and we require to rescale y ± so as to gain an airfoil of unit chord length. A rather coarse 64 × 16 C-type background mesh which extends about 20 chord length away from the airfoil is adopted for computation, see Fig. 9 for an impression of the mesh. The subsonic flow has configuration of Mach number M ∞ = 0.5, angle of attack α = 0 • and Reynolds number Re ∞ = 5000. We apply an adiabatic no-slip wall boundary condition along the airfoil, subsonic inflow and outflow boundary conditions at respective farfield boundaries. We employ also a combination of indicators η gi and η di with parameters ω 1 = 1.2, ω 2 = 0.3, ω 3 = 1.2, and ω 4 = 0.3. P 2 element as well as the first order IMEX-RK time discretization is used since this is a steady state simulation. The adaptive mesh and the computed Mach number contours are presented in Fig. 10, clearly showing that the elements along and behind the wake of the airfoil can be locally refined to capture the more complex flow field. The distribution of pressure coefficient along the airfoil wall is further displayed in Fig. 11 and is consistent with the results presented in [6,10].

Example 7 Supersonic flow past a NACA0012 airfoil
In this test case, the profile of the airfoil and the initial background mesh remain the same as in the previous subsonic case. The supersonic flow with Mach number M ∞ = 2 passes through the computational domain at an angle of attack α = 10 • and Reynolds number Re ∞ = 106. An adiabatic no-slip wall boundary condition is imposed on the airfoil and at the farfield part of the computational domain, supersonic/subsonic inflow and supersonic/subsonic outflow boundary conditions are assumed, respectively, based on the uniform freestream and the unit normal vectors to the outer boundary. A salient feature of this case is the formation of a bow shock in front of the profile. We will also apply a combination of indicators η gi and η di with parameters ω 1 = 1.2, ω 2 = 0.3, ω 3 = 1.2 and ω 4 = 0.3. Again owing to steady state simulation, P 2 element as well as the first order IMEX-RK time discretization is employed. Figure 12 shows the adaptive mesh, Mach number and pressure contours of the numerical solution. Evidently, the bow shock can be captured with a nice resolution by virtue of local grid refinement. The pressure distribution is also given in Fig. 13.

Example 8 Rayleigh-Taylor instability
In order to show the capability of our methods to capture the small-scale features of complicated flow structures, we consider the Rayleigh-Taylor instability phenomenon,   respectively. The source term ρ is added to the right hand side of the second momentum equation and ρv is added to the right hand side of the energy equation in (1). All the simulations are based on a 30 × 120 background mesh till final time T = 1.95, and the (k + 1)-th order IMEX-RK time discretization will be employed when P k element is used. For the NS simulations, we adopt a single indicator η gi with parameters ω 1 = 1.5, ω 2 = 0.5. Computational results with Reynolds number 50000 and 150000 are shown in Figs. 14 and 15, respectively. As can be seen, the interface between two fluids with different densities becomes thinner as the Reynolds number increases. Our adaptive methods can capture the interface and small-scale features precisely, and the resolution of P 2 element can even be comparable to that obtained with the ninth-order WENO scheme on a 600 × 2400 uniform mesh in [36]. The adaptive meshes and density contours of the Euler simulations with TVB limiter [3] and without TVB limiter are also presented in Figs. 16 and 17 by using adaptive parameters ω 1 = 1.2, ω 2 = 0.3. Due to the absence of physical viscosities, there exist more small-scale flow structures in these simulations, and the resolution of P 2 element can be comparable to that obtained with the ninth-order WENO scheme on a 480 × 1920 uniform mesh in [35].

Conclusions
In this paper, we have presented an h-adaptive LDG method for the solution of the unsteady compressible NS equations. To achieve the uniform high order accuracy both in space and time, avoid the extremely small time step restriction of explicit methods as  well as the construction of large Jacobian matrix in implicit methods, two semi-implicit time marching methods including the classical IMEX-RK methods and the SDC methods were also employed for temporal discretizations. In the process of mesh adaptation, indicators are computed in each element of the mesh based on the gradient of density, the divergence and curl of velocity field to pick candidate elements for refinement and coarsening. We also adopt the strategy of "refinement must, coarsening can" to determine the ultimate element lists for refinement and coarsening. These proposed methods have been Fig. 17 Euler simulations for Rayleigh-Taylor instability without TVB limiter. From left to right: adaptive mesh for P 1 element; 30 equally spaced density contours from 0.6548 to 2.4582 for P 1 element; adaptive mesh for P 2 element; 30 equally spaced density contours from 0.6852 to 2.3728 for P 2 element successfully implemented in a sequence of numerical examples, illustrating their expected optimal order of convergence, efficiency and capabilities in flow problems. Due to the high order accuracy and easy implementation, the h-adaptive LDG methods coupled with IMEX-RK or SDC methods provide an appealing alternative for unsteady compressible NS equations which require to accurately capture the detailed flow structures at any time.