Shock capturing with the high‐order flux reconstruction method on adaptive meshes based on p4est

High order schemes have been investigated for quite a long time, and the flux reconstruction (FR) scheme proposed by Huynh recently attracts the attention of researchers due to its simplicity and efficiency. Building the framework that bridges discontinuous Galerkin (DG) and spectral difference (SD) schemes, FR recovers DG and SD conveniently with a careful selection of parameters. In this article, FR scheme is realized based on the framework of p4est, an open source adaptive mesh refinement library. The shock capturing ability of localized Laplacian artificial viscosity and in‐cell piecewise integrated solution methods are compared. Curved boundary treatment for high order schemes is adopted. The performance of developed code is estimated in both one and two dimensions including curved boundary and shock cases, and some attractive results are obtained.

and viscous cases [11][12][13] on structured and unstructured grids. 14,15 According to the property of different parameters, some optimization works have been done. 16 Also, FR scheme has been applied in real applications like large-eddy simulation (LES), 17 overset grid, 18 multi-phase simulations, 19,20 and so forth.
It is convenient for high order schemes to maintain high accuracy in smooth flows. However, capturing of discontinuities remains challenging in DG type methods since non-physical oscillations could occur near these regions, which may cause fake solutions or lead to divergence. Among methods suppressing non-physical oscillations, the limiter and artificial viscosity are well-known techniques. Limiter tends to limit solution to a specified range which acts directly to conservative variables rather than flux. In DG type schemes, Cockburn and Shu 2 proposed a TVB type limiter, Du et al. 21 developed a WENO type limiter for CPR scheme. Park et al. 22 introduced a multidimensional limiter for CPR. Li and Wang 23 developed accuracy preserving limiter. Zhang and Shu 24 proposed the maximum-principle-satisfying limiting procedure for high order schemes. The localized artificial diffusivity (LAD) 25 and localized Laplacian artificial viscosity (LLAV) 26 methods are more commonly adopted in artificial viscosity methods. Moreover, there are other shock capturing methods like filtering technique 27 and the newly developed in-cell piecewise integrated solution (IPIS) method. 28,29 Similar works in DG scheme has been done with a DG/FV style to perform shock capturing by Sonntag and Munz 30 and Dumbser et al. 31,32 There are schemes in different framework like the arbitrary Lagrangian-Eulerian 33 (ALE) or invariant domains method 34 (IDM) that provide different perspectives. It is of great significance to clarify the differences and characteristics of different shock capturing methods, which can guide the follow-up work. In this article, we test and compare the performance of LLAV and IPIS methods.
Researchers have realized that in DG type schemes false solutions may occur 35,36 when curved boundary exists in wall geometry as the flux points may not locate exactly on true boundary. In this article, the Coons curve reconstruction adopted by Hindenlang et al. 37 is applied to perform boundary fitting. The research of Coons curve reconstruction could be traced back to the research of Gordon and Thiel. 38 And recent work of Heltai et al. 39 builds a framework that adds new nodes with only two operations. In current research we do not focus on the improvement of the algorithm, and only a simplified application of boundary treatment is adopted, more detailed information could be found in referenced researches.
The adaptive mesh refinement (AMR) technique has been widely adopted to improve computational efficiency with much less computational cost compared with uniform global refinement to the same level with AMR. AMR adjusts grids according to the defined criteria concerning different flow properties. Generally, three types of AMR 40,41 exist currently, the adaptive h-, r-, and p-refinement methods. P4est, 42 the scalable h-adaption AMR library adopted in our research, eliminating the need for developers to consider complex mesh adaption algorithms.
In this article, we combine the advantages of both FR scheme and AMR technique to improve solving performance. The LLAV and IPIS methods shock capturing methods are realized and their performances are compared. The remainder of this article is organized as follows. Section 2 introduces FR scheme, curved boundary treatment, and time integration method. Different shock capturing methods and AMR technique are introduced in Section 3. Section 4 presents and discusses the numerical results obtained. The efficacy of shock capturing methods is compared in both one and two dimensions. The accuracy and efficiency of AMR are also evaluated in this part. The research is summarized in Section 5.

One-dimensional FR scheme
In this section, FR scheme is introduced in one dimensional conservation law First consider partitioning the computational domain [a, b] into N non-overlapping cells, and define ith cell with To normalize FR scheme, each cell is transformed into a standard cell [−1, 1], and the local coordinate is , the transform relation between x and is ∕2 is the coordinate of the center of ith cell. K points k , k = 1, 2, … , K are selected as solution points to define a standard cell, the points on the left and right boundaries are flux points. Assume the value u i,k (t n ) at time t of each solution point is known, without confusion, write u i,k (t n ) as u i,k (t). In FR, there are different solution points can be used, such as Gauss points, Lobatto points, or equidistant points. With Lagrange interpolation, a polynomial of degree K − 1 can be constructed as follows.
where k is the kth Lagrange basis function The flux distribution on the boundary of adjacent cells is generally different. To construct a continuous flux function, Huynh introduced the correction concept, in which a correction function g( ) of degree K is added to flux function.
To consider the interaction of neighboring cells, common values on the boundaries need to be computed. Use L and R to represent the variables on left and right boundaries, here the flux function is not continuous. The Roe Riemann flux is adopted as common flux.
With correction function g( ), continuous flux function can be written as here the last two terms on the right side of the equation are the correction term on the left and right boundaries. g L and g R are the left and right correction functions, respectively, which satisfy There are different correction functions can be selected, which can be found in Huynh's research. Radau function makes the scheme be equivalent to DG scheme. In this article, Radau function is adopted to manage high accuracy and according to Huynh, 6 this leads to accuracy of order 2K − 1 in linear systems for long-term numerical simulations.
The following differentiated flux function is obtained next to perform time integration.
By now, the derivative of flux in the computational domain is obtained. With the transform equation of x and , the derivative in the physical domain could be obtained And the conservation law becomes Huynh pointed out that in linear cases the choice of correction function has a significant influence on stability and precision, the order of the scheme is irrelevant with the type of solution points. Jameson et al.'s 43 research, however, indicates that the type of solution points is important and can affect the energy property and aliasing error when the conservation laws are nonlinear. Gauss points can eliminate the aliasing error and are chosen as solution points in current research. Explicit or implicit time marching methods can be applied to solve the differential Equations.

Two-dimensional FR scheme
For a mesh constructed with quadrilateral cells, FR scheme can be extended to two dimensions without much effort when the grid lines of applied mesh are all horizontal or vertical. The calculation can be conducted in a dimension by dimension manner.
For arbitrarily shaped cells, the coordinate transformation between physical and computational domains is necessary. The transformation from a quadrilateral cell to a standard cell is shown in Figure 1 and can be expressed as where M i ( , ) is the shape function of each node and N is the number of nodes that define the standard cell, x i , y i is the global coordinate of each node. In this article, two-dimensional Euler equations are adopted. And the conservation laws are where |J| represents the determinant of Jacobian matrix.

Curved boundary treatment
For high order schemes, it has been investigated and known that curved boundary needs to be treated with additional effort. 35,36 In FR, if the boundary is represented with straight lines, the flux points at wall boundary may not represent the geometry accurately. This could lead to false solutions or divergence issues.
Only several software supports high order mesh generation yet, curved boundary treatment often needs to be accomplished by researchers themselves. In general, four-node standard cell is suitable for cells constructed with straight lines, while the six-node cell or eight-node cell shown in Figure 2 should be considered when curved boundary exists. The six-node cell has a boundary edge of third order which is sufficient for 2D wall fitting. In viscous cases, the cells near the boundary layer are highly anisotropic and the first layer of cells is very thin in wall-normal direction. The curved boundary could intersect with one or more layers of boundary cells if only six-node cell is applied. This could be avoided using the eight-node or other curve-edged standard cells. In current research, only inviscid tests are performed, and the six-node cell would meet the demand.
In real applications, it is always difficult to obtain the exact shape function of wall boundary, thus the high order function constructed with local information is often preferred. In this article, local Coons curve reconstruction 37 is adopted. Figure 3 illustrates an example of the reconstruction procedure, where A, B, and C are three linear cells, the normal vector and the length of each boundary face are n A , n B , n C , and l A , l B , l C , respectively. n 1 and n 2 are normal vectors of vertex 1 and 2 computed with weighted averaging. For example, Not all normal vectors of boundary nodes need to be computed with two edges sharing the node. Averaging would cause over-fitting when the angle between two edges is sharp, such as the trailing edge of an airfoil, and the normal vector remains unchanged under this circumstance.
Coons curve reconstruction is presented as follows. Take cell B for example, with the normal vectors of vertex 1 and 2, third order Hermite interpolation is applied to construct the boundary function.
F I G U R E 2 Cells with straight edges and curved edges

F I G U R E 3 Curve fitting diagram
The Coons curve fitting only involves the information of neighboring cells, which makes coding easier and has some advantages in parallel computing. In static grid system, the curved boundary treatment needs to be performed only once, and the time saved is limited compared with simulation time. In cases with deformable geometry, the procedure should be conducted along with each deformation, and this would reduce some time in boundary treatment.
For a four-node standard cell, the iso-parametric mapping function is And the mapping of the adopted six-node standard cell is

Time integration method
The three stage strong stability preserving Runge-Kutta method 44 is adopted for time marching. When high order spatial discretization is coupled with high order time discretization, it preserves the stability of simulation and has low memory consumption. For conservation laws The Runge-Kutta method is written as here h is the cell size, is the characteristic speed.

SHOCK CAPTURING AND AMR
Shock capturing methods are applied to suppress the Gibbs phenomenon generated near discontinuities. In this section, we introduce and evaluate the performances of the artificial viscosity and IPIS methods. The AMR technique is adopted to improve the computing efficiency.

Artificial viscosity
Artificial viscosity technique has been studied for a long time in high order schemes. But how to define the amount of viscosity is still an open question since there are adjustable parameters that are problem dependent. LLAV method proposed by Persson and Peraire 26 has the property of limiting shock within a cell and was modified by Yu et al. 45 In this article, the modified version is adopted. The LLAV method adds a Laplacian diffusivity term in control equations as shown below On the right-hand side is the Laplacian artificial viscosity term, is the viscosity parameter, which is dependent to the property of local flow. In the research of Persson and Peraire, ∼h∕p, h is the cell size, p is the order of solution polynomial. And can be expressed as follows here s e is the smoothness indicator which will be introduced later. In Equation (22), 0 is the artificial viscosity which can be regarded as the function of Δ max and Péclet number Pe, size of cell h and the maximum characteristic speed | | max here Δ max ∈ [0, 1] is the ratio of the distance between two solution points and cell length. There are different selections of functions for f (Δ max , Pe). 45 The linear function is used in current research.
The viscous term in the equation is treated with the following method. 26 With the auxiliary variable q, the system equations can be treated similarly as traditional FR scheme. The common flux of viscous terms is computed by averaging the flux of neighboring cells. The amount of is set constant inside cell, although the constant artificial viscosity may lead to some unstable issues, but in current research the results show that constant is acceptable with smaller CFL number. And in the implement of AMR technique, with the volume of cells in shock region getting smaller, the spurious phenomenon generated by non-smooth tends to be weaker.
Persson introduced the shock indicator below to help remove artificial viscosity in smooth regions. First, the solution polynomial is expressed with orthogonal basis Here choose density as u. N(p) is the number of basis functions, i is the orthogonal basis, u i is the coefficient for each basis function. And for each element, the truncated P − 1th order polynomial iŝ And the smooth indicator is where (*,*) indicates inner product in L 2 and the parameter S e ∼1∕p 4 , s 0 = −3 log 10 (p), in Equation (22) is an experimental coefficient. Based on the research of Yu et al., 45 has a larger influence compared to other parameters, and small value of could not limit spurious oscillations effectively while the large value smears the discontinuities excessively. And should be set between 3 and 6.

In-cell piecewise integrated solution
The IPIS method is the shock capturing method proposed by Huera et al. 28 in DG scheme and was extended to FR by Yi et al. 29 The idea of IPIS is to construct constant solution in each sub-cell and reconstruct the distribution function. It does not need information from neighboring cells compared with artificial viscosity method. To construct the piecewise constant solution, each cell is partitioned into sub-cells containing solution points. For FR scheme, solution points locate inside the cell, sub-cells could be constructed containing each point as shown in Figure 4. For 2D cases, the number of sub-cells is the same with the number of solution points. If 3 × 3 solution points are in each cell, then 9 sub-cells containing each solution point respectively will be defined. The constant solution at each solution point is constructed as follows

F I G U R E 4 Sub-cell construction in 2D
here e j is the sub-cell containing jth solution point, V j is the volume of e j . New solution u new e,j at solution point j is obtained as follows where ∈ [0, 1], which depends on the smoothness of flow field. And it can be calculated the same with the artificial viscosity coefficient in LLAV.

AMR method
There are finite element packages support AMR like deal.II, 46 MFEM, 47 DUNE, 48 Firedrake, 49 and so forth, which provide useful functions to assist users to perform researches. Current code is developed under the framework of p4est to realize AMR efficiently. P4est 42 is an open source AMR library which provides abundant API for different users. It has high parallel efficiency and is scalable. The users need to define the adaption criteria and data projection while grid is refined or coarsened. We have adopted P4est in our former works in which AMR works well, moreover we are familiar with the data structure of it. We would like to develop a FR solver that focus on the realization. This would accumulate experience on the connection of different processes in details for further research.

Grid adaption criteria
In AMR, in addition to the complicated data structure and algorithms, another key aspect is the determination of where to refine or coarsen cells. In fluid dynamics the commonly adopted criteria is based on the local gradient of physical variables like density, temperature, pressure, entropy, and so forth.
Different indicators should be used for different flows. In flows containing discontinuities, the shock detector adopted in shock capturing could be used as indicator. When s e < s 0 − , which means solution is smooth and the cell is marked "coarsen" unless the minimum user defined level is reached. Otherwise it indicates that discontinuity exists, and the grid is marked "refine." In flows without discontinuity, vortex property is more important. The following criteria presented in Reference 50 are used for vortex flows.
where V is the speed vector, N is the number of cells. If ci > c , mark the cell "refine," otherwise it is marked "coarsen." In current work, the refine and coarsen processes are carried out one after another. After refinement process, the indicators of child cells are set to a value larger than refine threshold instead of projecting. Coarsen step is conducted after refinement process. Each cell will go through both refine and coarsen once respectively in each adaption process, which is not a recursive procedure.

3.3.2
Nonconforming data projection FR scheme needs to calculate the common flux on cell boundaries at each time step. The common numerical flux can be calculated directly when computational grid is uniform and stationary as the edges and flux points coincide between neighboring cells. When AMR is conducted, the common edges may not share single edge, called non-conforming edges. Under this case, the flux points on each edge are not collocated, as shown in Figure 5. To calculate common flux, the mortar element method (MEM) is applied, a more detailed description can be found in References 51 and 52. And L2 projection is used to project data when mesh is refined or coarsened.

NUMERICAL RESULTS
To test the shock capturing ability of artificial viscosity in LLAV and IPIS, we have tested shock cases in both one and two dimensions in Section 4.1. The shock capturing effect with AMR is tested with different cases, which is shown in Section 4.2. The time marching method adopted is three-stage strong stability preserving Runge-Kutta method introduced in Section 2.4.

One-dimensional SOD shock tube
The SOD shock tube problem is a benchmark case for shock capturing. It is initialized with a Riemann condition and three different flow phenomenon namely the rarefaction wave, contact discontinuity, and shock wave could be generated and it has been widely used in the examination of numerical schemes. The initial condition is The grid size is dx = 1/100. Figure 6 represents the results with artificial viscosity and IPIS methods, CFL = 0.5. They are compared with reconstruction degree of 2 and 3 (third order and fifth order) at t = 0.2 s. The results reveal that both methods have suppression effect on spurious oscillations. Both methods capture shock with good accuracy while the dissipation to contact discontinuity is slightly larger. IPIS shows smaller dissipation, which makes both the contact discontinuity and shock profiles shown in Figure 6 sharper compared with LLAV. In addition, as the solution order increases, the numerical solution gets more precise, which is consistent with expectation.
Additionally, to demonstrate the effect and the proper range of the parameter , we did tests from = 2 to 10 using LLAV with fifth order accuracy. As Figure 7 depicts, controls the range where shock capturing methods get active. It is obvious that when is small, the shock becomes highly oscillatory while the large ones smear the shock. From the result, the range of ∈ [3,6] is suitable for current research and this is in agreement with the research of Yu et al. 45

4.1.2
One-dimensional Shu-Osher problem Shu-Osher problem solves the interaction between a flow with sinusoidal density distribution and a shock moving with speed of Mach 3 to the right. It contains both small flow structures and strong discontinuity which is challenging for high order schemes. The initial condition is The computational domain is [−5, 5], with grid size of 1/20, and the order of accuracy is 3, Figure 8 presents the density distribution at t = 1.8 s. The parameter is set 4, CFL = 0.2. The reference result is calculated with fifth order WENO scheme using 2000 cells. The results indicate that both methods solve the problem well and there is little oscillation near shock. It could also be observed that IPIS still shows smaller dissipation and dispersion than LLAV.

Supersonic wedge
The supersonic flow over a wedge is used to test the performance of methods in 2D, the geometry and grid are displayed in Figure 9, mesh contains 95 × 31 cells, the order of accuracy is 3 and the initial condition is ( , u, v, p) =  Figures 10 and 11 give the cloud pictures and isolines with two shock capturing methods.
The results indicate that LLAV solves the Mach stem at the upper boundary better. LLAV shows more dissipation in 1D and is compatible in 2D compared with IPIS. In addition, LLAV runs slower but more robust during simulations. In current tests, IPIS spends about 30% time lower than LLAV. And the memory consumption of LLAV is about 10% higher than IPIS. The reason is that the LLAV needs to store the viscous terms with additional space. And the treatment of viscous term is more time consuming than the solution reconstruction acts directly to conservative variables, but this might be an advantage in viscous flows as the space of viscous terms could be used directly or physical viscosity could be used as artificial viscosity and this still needs to be tested.

Cases with AMR
The cases in this section are tested with LLAV shock capturing method and degree 2 reconstruction (third order) is adopted. For all tested cases, the maximum level of refinement of AMR is two.

Double Mach reflection
The Double Mach reflection problem 53 is used to test the ability of schemes to capture strong shock. The flow is initialized by one shock with the speed of Mach 10 reflected by a solid wall with incident angle of 60 • . The computational domain is [0, 3] × [0, 1], the lower boundary at 0 ≤ x < 1∕6 and the left boundary are set post-shock initial condition, the lower boundary at 1∕6 ≤ x ≤ 4 is set slip wall boundary condition, the right boundary is set free outflow condition, upper boundary is set the same as the exact movement of shock. is set 5, CFL = 0.1. The pre-shock and post-shock condition is In this case, initial grid size is 1/30, order of accuracy is 3, and artificial viscosity is applied for shock capturing. Figure 12 shows the density contour and mesh distribution at t = 0.2 s. Figure 13 shows the zoomed in results of double Mach reflection. The results show that artificial viscosity method is capable of capturing strong shock. The grid is refined in shock region and vortex flow structure in Mach stem region can also be obtained.

Transonic NACA0012 airfoil
In this part, the transonic flow passing NACA0012 airfoil is tested with inflow condition Mach = 0.8, angle of attack = 1.25 • . This test has been widely adopted as the benchmark case in testing the accuracy of numerical schemes. Initial computational grid contains 2400 cells, and curved boundary treatment is applied. is set 4, CFL = 0.4. Figures 14 and 15 display the mesh distribution and results without mesh refinement. Figures 16 and 17 show the results with mesh refinement and final mesh distribution. Pressure distribution is compared with the results of Jameson. 54 and Yang et al., 55 which is presented in Figure 18. It is clear that the shock profile gets sharp with AMR and the results are on good agreement with their works. Mach contour reveals that shock on upper surface becomes sharper with mesh refinement and the range of artificial viscosity distribution is narrower. Based on pressure coefficient comparison, when shock is weak, like the shock at lower surface, the mesh refinement plays a smaller role. And the overall cell number maintains approximately 2700 during the calculation process. Although mesh refinement could lead to a smaller time step, but the computation time compared with global refinement to obtain compatible accuracy will be much less. This will be discussed in next subsection.

Efficiency verification of AMR
We compared the computational efficiency of AMR versus uniform mesh refinement. In this case, the maximum refine level is 2 for AMR part, and all cells of initial mesh are refined to level two to form globally refined mesh. During simulation number of cells with AMR remains around 2700 while the globally refined mesh contains 38,400 cells. The results are compared with the result of Yang et al. 55 The simulation stops when |lift − YangLift| <= 0.005. The time and lift/drag comparison are displayed in Figure 19. It is clear that the time consumed is extremely long with globally refined mesh and the simulation with AMR only takes 10% of the time spent in globally refine mesh. The supersonic wedge is tested again to demonstrate the mesh adaption strategy, simulation stops at time = 3 s. In this part, the mesh is tested with the maximum refine level from 0 (referring to initial mesh) to 2. Figure 20 depicts the mesh variation under maximum refine level 2. It is clear that the cells in shock region are refined, the refine and coarsen area changes along with shock wave motion. In Figure 21, the results from initial mesh to 2 level refinement are presented, the shock gets sharper and the smearing effect is alleviated. In current tests, no obvious spurious oscillation was observed with AMR.

CONCLUSION
In this article, the FR scheme is realized in one and two dimensions under the framework of the AMR library p4est. The shock capturing capabilities of the LLAV and IPIS methods are evaluated and compared. Both methods are suitable for shock capturing. IPIS method is relatively simpler and has smaller dissipation in one-dimensional. Both methods resolve shock wave better compared with contact discontinuity. Meanwhile, LLAV takes about 30% longer time and 10% larger memory consumption, but it works more robustly. Coons curve reconstruction is applied to fit curved boundaries. The test of NACA0012 airfoil confirms the correctness of the method. With the help of p4est, AMR is integrated into our code. The tested cases reveal that the results with AMR capture more details than those without. The adaption region varies with flow properties. The smearing effect of shock capturing methods is alleviated with AMR and no obvious spurious effect is observed during simulations. It can be concluded that AMR technique is efficient in improving solution efficiency.
One disadvantage in AMR with explicit time marching method is that the physical time step would get smaller with minimum cell size decreases. This would lead to a longer time consumption according to the maximum refine level. Thus the accelerate technique or the implicit time marching method will be the direction of our further research.