A characteristic-featured shock wave indicator for simulating high-speed inviscid flows on 3D unstructured meshes

In this work, we extend the characteristic-featured shock wave indicator based on artificial neuron training to 3D high-speed flow simulation on unstructured meshes. The extension is achieved through dimension splitting. This leads to that the proposed indicator is capable of identifying regions of flow compression in any direction. With this capability, the indicator is further developed to combine with h-adaptivity of mesh refinement to improve resolution with less computational costs. The present indicator provides an attractive alternative for constructing high-resolution, high-efficiency shock-processing methods to simulate high-speed inviscid flows.


Introduction
Recently, high-order, high-resolution schemes are widely used in computational fluid dynamics for their strength in capturing detailed features of wave structures. Specially, Discontinuous Galerkin (DG) methods [1][2][3][4] can treat wave propagation well due to its upwind numerical flux which solves a Riemann problem arising from the discontinuous representation of the solution at element interfaces. As compared to finite difference (FD) and finite volume (FV) methods, DG schemes are easier to achieve high order and easier to handle complex boundaries on unstructured meshes. However, DG methods have a higher probability of making solutions suffer from unphysical oscillations as compared with low-order methods. The oscillations are produced by over-fitting shock wave discontinuities using high-order degrees of freedom (DoFs), and usually cause robustness issues of simulating high-speed flows. In addition, DG is recognized as more expensive in terms of computational costs for that it requires more DoFs to achieve the same accuracy as compared to classical finite element (FE) methods.
To eliminate unphysical oscillations near shock/contact wave discontinuities and improve robustness of DG methods on simulations of high-speed flows, a number of shock-processing (correction) techniques are designed to correct high-moment DoFs between spatial discretization and temporal discretization, such as limiting techniques [3,5], solution reconstruction [6][7][8], and artificial diffusion [9][10][11]. However, the use of correction techniques usually leads to inefficiencies of numerical schemes or even deterioration of accuracy. Therefore, it is necessary to employ a troubled-cell or shock wave indicator to guide accurate application of correction techniques. There are several practical troubled-cell indicators, such as the classical Total Variation Bounded (TVB) based indicator [3], KXRCF shock detector [12], Harten's troubled-cell indicator [13], the posteriori subcell-based indicator in [14], and a simple indicator in [15]. There are still very recent Artificial Neural Network (ANN) based troubled-cell indicators in [16][17][18] to improve the accuracy of discontinuity capturing. However, these indicators mentioned above mainly detect troubled-cells based on the posteriori errors of numerical oscillations, usually cannot theoretically avoid detecting noise such as wrongly marked extremums, and part of them need increased computational costs to evaluate their values.
In our recent work, we employ an Artificial Neuron (AN) with characteristiccompression-imbedded data to propose a characteristic-featured shock wave indicator. The resultant indicator is concise, compact and explicable, and we have proved in [19] that its output can measure the compression of characteristic curves, which evolve into shock wave discontinuities. The indicator has been extended to 2D unstructured meshes and its performance combined with various correction techniques was presented in [20]. In this work, we develop this indicator to 3D unstructured meshes via dimension by dimension. This extension technique makes the present indicator capable of detecting compressible waves in any direction. Such a feature enables the present indicator to work closely with hadaptivity of mesh refinement to further improve resolution and efficiency of simulation of high-speed inviscid flows.
The rest of the paper is arranged as follows. In Section 2, we briefly introduce the governing equations, spatial Discontinuous Galerkin (DG) discretization, temporal discretization and correction techniques including artificial viscosity. The description of the characteristic-featured shock wave indicator, development on 3D unstructured meshes, and some corresponding properties are presented in Section 3. Mesh refinement strategy based on the present shock wave indicator is presented in Section 4. Numerical experiments are provided in Section 5 to show the performance of the present indicator combined with various correction techniques and mesh refinement on multi-dimensional unstructured meshes. We make a few conclusions in Section 6.

Governing equations
We consider a multi-dimensional conservation law where U ∈ R m is the vector of conserved variable and F(U) = (F d (U)) ∈ R m × R d . m, d represent the dimension of conserved vector and spatial space, respectively. In this work, we focus on the Euler equations, where the conservative variable vector U and the advection flux vector F(U) in 3D space are here ρ, p, E denote the density, pressure, and total energy of the fluid, respectively. u = (u 1 , u 2 , u 3 ) is the velocity vector of the flow in the coordinate direction. The pressure can be computed from the equation of state Here, γ = C p /C V denotes the ratio of specific heats. We set γ = 1.4 for all the test cases in this work.

General DG formulation
The Discontinuous Galerkin (DG) scheme is used to discretize the weak form of the governing Eq. (1) over the domain where (= ∂ ) denotes the boundary of and n d is the unit outward normal vector to the interface boundary. We assume that the domain is subdivided into a collection of non-overlapping elements e . We introduce the following broken space of vector-values polynomials with a degree of p Provided the basis functions {φ e i (x)} of V p h , the local numerical solution U h on e can be computed by the basis function The Taylor basis functions [21] are used in this work. For sake of simplicity, let us consider a quadratic polynomial solution U h in 3D space which is expressed as with 2 2 y 2 d , Replacing test function φ in (4) by elemental basis function φ e i (x), we can obtain the semi-discrete form by applying weak formulation on each element e , where i is the index of polynomial basis on each element e , F d is a consistent numerical flux.

Temporal discretization
The spatial discretization (9) leads to a system of ordinary differential equations, where M denotes the mass matrix, U is the solution vector, and R(U) is the residual vector with the form of A high-order and stable explicit temporal discretization SSP-RK3 [22] scheme is applied to treat unsteady problems. This work focuses more on steady problems, in order to speed up the convergence and improve the stability of simulation for steady problems, implicit temporal discretization is presented in this work. Using Euler implicit time-integration, the spatial discretized Eq. (9) can be linearized in the temporal direction and written as where ∂R ∂U is the Jacobian matrix, t is the time increment, and U n = U n+1 − U n is the solution difference between time level n and n + 1. This system of linear equations at each time step is solved by the GMRES method with a LU-SGS preconditioner [23].

Correction techniques
As mentioned in the introduction, there are a number of correction techniques to treat shock wave discontinuities. For example, we employ the Barth-Jespersen limiter [24] in an explicit temporal scheme to treat unsteady problems in [20], and employ the artificial viscosity with strong-residual term in an implicit temporal scheme to treat steady problems. Here, we describe artificial viscosity in details.
The artificial viscosity technique can treat shock waves through adding an adaptive artificial viscosity term to bring parabolic properties into hyperbolic conservation laws (1) to smear the shock wave discontinuities. Thus, the semi-DG scheme is modified into Parameter e in this work is a Hartmann-type coefficient of strong residual where C and β are problem and mesh dependent parameters which are usually set to 0.0001 − 0.01, 0.01 − 0.5, respectively.

The characteristic-featured shock wave indicator
In our recent work [19,20], we proposed a characteristic-featured indicator based on training an artificial neuron (AN) to detect shock wave discontinuities. We firstly describe the indicator on 1D mesh, next develop the shock wave indicator on 3D unstructured meshes, then show some properties of the present indicator, and finally provide its application prospects as a shock wave indicator.

Indicator on 1D mesh
For 1D scalar conservation laws ∂u ∂t we consider a numerical solution u h (x, t) on a perturbed mesh, h e denotes the meshsize of element I e ,ū e = 1 h e I e u h (x, t)dx is the cell-avaerge of variable u inside I e . The characteristic-featured shock wave indicator through training an artificial neuron (AN) [20] using the information of neighbours (I e−1 , I e , I e+1 ) is presented as follows, Here, W = 10, M 1 = 4, M 2 = 1, andλ L ,λ R are the left and right side weighted averages of eigenvalue λ(u) = f (u) which are convex combinations of cell-averagesλ e−1 ,λ e and λ e+1 ,λ e , respectively (Refer to Fig. 1), in [20],λ L ,λ R are set to integral-averages In practice, we found M 1 in (15) is not sensitive to the convergence of AN. We retrain an AN with 3 argumentsλ L ,λ R , h M (h M := max {h e−1 , h e , h e+1 }) and without hidden layer, and the AN convergents as well and maintains a similar working mechanism. In order to make the indicator more concise, we adopt a model with less input arguments and keep most coefficients in (15) unchanged. The expression of the modified shock wave indicator is present as follows, Here, W = 10, M 2 = 1 keep unchanged, M 1 = 12,λ L ,λ R are defined by (16) and h M = max {h e−1 , h e , h e+1 } is the local maximum mesh-size. We drop the subscript "M" in the following text for convenience.
In [19], we have proved that the output of indicator (17) measures the admissible jump values of eigenvalues, the greater the value ofλ L −λ R is, the stronger the admissible jump is (a big admissible jump value implies characteristic compressing and leads to the development of a shock wave), at the same time the output of the indicator (17) is closer to 1.

Extension to system
For 1D system of conservation laws (1), the indicator is extended by characteristic field splitting, that is to apply the scalar indicator (17) to each genuinely nonlinear characteristic field to detect the compressing or intersection of characteristic curves (refer to [25]). The expression is presented as follows, Here,λ i,L ,λ i,R are the left and right side-weighted cell-averages of eigenvalue λ i of i-th genuinely nonlinear characteristic field, respectively, and out e (λ i,L ,λ i,R ; h) is defined in (17). Especially, for the 1D Euler equations, the indicator (18) becomes where a = γ p ρ represents the speed of sound. It should be noted that the contact discontinuity can be detected with the indicator (19), because there is always an admissible jump occurring to the 1-or 3-eigenvalues across the 2-contact discontinuity.

Development on 3D unstructured meshes
The indicator (17) is extended to 3D unstructured grids through dimension splitting, that is to construct side-weighted averages of eigenvalueλ on each dimension, and then apply (17) in each spatial dimension to accomplish the indication. Note that compression of characteristic cones in any direction results in compression of the characteristic curves projected onto the x or y or z-direction. Therefore, multi-dimensional indicator through dimension splitting can capture shock discontinuities aligned to any directions, although it is only applied along x, y and z-directions.
Thus, careful evaluation of the weights of neighbours to constructλ becomes crucial to avoid the influence of mesh shape and mesh size. Here, we present a relatively simple construction method shown in our recent work [20].
We set the 3D Euler equations on the following mesh (refer to Fig. 2) as an example to show the extension procedure. Firstly, we present the following notation to be used, 1 , 2 , 3 are the neighbours of 0 , b k = (x k , y k , z k ) is the centroid of k , | k | represents the volume of k , A k represents the volume ratio with the form of | k | | 0 |+| k | , n k = (n kx , n ky , n kz ) represents the unit vector of centroid connection vector b 0 b k |b 0 b k | , λ k represents cell-averages of eigenvalues in k in the x, y, z direction, λ = u 1 ±a, u 2 ±a, u 3 ± a.
Based on the methods shown in [20], in order to construct side-weighted averages of eigenvalueλ and reduce the influence of mesh shape and mesh size, we need to solve two key issues, (a) one is how to select the left or right side element of 0 on each spatial dimension among all neighbours, such as right-side element of 0 in the z-direction (named as top-side element); (b) the other is how to select the "weight" for each sideelement, such as how to evaluate the weights for each top-side element when constructing weighted averageλ T . The procedure of constructing side-weighted averagesλ on each spatial dimension is presented as follows.

Fig. 2 Geometry information for 3D Euler indicating on unstructured grids
Step-1: Selection of left, right side-element on each dimension. The direction of n k is used to determine the position of k to 0 . k is recorded as the backward element of 0 if n kx < 0, similarly, For example, refer to Fig. 2, 1 is marked as a forward element in the x-direction, a left-side element in the y-direction, and a top-side element in the z-direction.
Step-2: Evaluation of weights of each side-element. As the position of k to 0 on each dimension is determined, the next is to evaluate the weights that λ k contributes for constructing side-weighted averagesλ on each dimension. The weights are products of two parts, the direction weight and the size weight, (a) direction weight. The direction weight represents the contribution of k on each spatial dimension, for n 2 kx + n 2 ky + n 2 kz = 1, the greater n 2 kx is, the more contribution k makes in the x-direction, while the less contribution k makes in the y, z direction. Therefore, n 2 kx is selected as the direction weight for constructing backward-weighted averagē λ b or forward-weighted averageλ f , and n 2 ky , n 2 kz are selected as the direction weight for constructingλ L orλ R ,λ B orλ T , respectively.
(b) size weight. The mesh size needs to be considered as another factor for weights evaluating. The greater the volume | k | is, the more contribution k makes on each dimension. Therefore, we set the volume ratio A k = | k | | 0 |+| k | as the size weight for constructingλ on each dimension.
Based on Step-2 and refer to Fig. 2, we can construct the forward-weighted averageλ f 1 contributed by 1 , which is a convex combination ofλ 1 andλ 0 , where n 2 1x is the direction weight and A 1 is the size weight for 1 in the x-direction, w f 1 is normalized weight combined by both direction and size weights.
Similarly, we havē At last, we assume all the same side-elements on each spatial dimension have an equivalent opportunity to contribute to the construction of side-weighted averages. For example, refer to Fig. 2, only 1 is determined as a left-side element in the y-direction, both 1 and 2 are determined as top-side elements in the z-direction, thus we havē In addition, we define the characteristic mesh-size h x , h y , h z to be used as follows, According to the steps above and combined with formulae (20), (21) and (22), sideweighted averages of eigenvalueλ b ,λ f ,λ L ,λ R ,λ B ,λ T can be constructed easily. Replacing h in (17) by h x in the x-direction, h y in the y-direction, h z in the z-direction, andλ L ,λ R in (17) byλ b ,λ f in the x-direction,λ L ,λ R in the y-direction, andλ B ,λ T in the z-direction, we can obtain a concise and generalized characteristic-based shock wave indicator on unstructured grids as follows, Here, out e is the indicator for system equations (18),λ b ,λ f ,λ L ,λ R ,λ B ,λ T are defined by (20), (21), (22). (24) is used as the shock wave indicator in this work. Based on statements in Section 3.1 , we can conclude that (24) can detect shock wave discontinuities caused by compressing or intersecting of characteristic. Especially it can detect contact waves in the Euler equations.

Indicating procedure based on the present indicator
We present two approaches to detect troubled-cells having shock wave discontinuities based on the present indicator (24) on arbitrary grids. One is to flag cells by a critical value s, that is out 3D e ≥ s =⇒ e is flagged as a troubled cell. out 3D e < s =⇒ e is not treated.

(25)
The parameter s is set to 0.5, while s can be set smaller for identifying weaker admissible jump strength. This approach is suitable to be combined with TVB or WENO limiting techniques when using explicit temporal schemes to simulate unsteady flows.
The other is to flag cells with a fixed fraction θ. Only the top θ fraction of elements with output of (24) is flagged. This approach is suitable to be applied with artificial viscosity in implicit temporal schemes to treat steady flows, and suitable to be applied in adaptive mesh-refinement strategies as well. This use will be described in the next section.

Mesh adaptive method
Adaptive methods provide an attractive alternative for reducing computational costs of DG schemes. Appropriate adaptive methods are able to focus computational resources on needed regions to further capture the detailed features of wave structures (such as shock and contact waves, vortex) and avoid the additional increase in computational costs caused by global mesh refinement or global order enrichment.
The essence of adaptive methods is how to flag elements which need focusing on. As mentioned in the last section, the indicator (24) is able to measure the compressing of characteristics. This feature results in the indicator being naturally suitable as a special error estimator which measures the unidirectional variation of eigenvalues to capture specific wave structure such as shock, contact wave and vortex through detecting the compressing of characteristics.
There are two strategies for mesh refinement once the local error estimator is determined. One is adaptively refining the mesh with a fixed fraction (each time only the top or bottom fixed percent of mesh is refined or coarsened); the other is refining the mesh with a fixed number (each time only a fixed number of mesh is refined or coarsened). In this work we select the former strategy based on the following error estimator for mesh (quadrilateral grid) refinement, Here, the subscript d denotes the dimension of space, out d e is the 1D indicator (17), λ L,d ,λ R,d are side-weighted averages of eigenvalues in the d-direction defined by (20), (21), (22), and h d is the local characteristic mesh-size in the d-direction defined in (23).
The whole adaptive algorithm based on the above error estimator (26) can be briefly summarized in Algorithm 1 as follows: The performance of the proposed adaptive method is presented and discussed in Section 5.

Numerical results
We now provide the performance of the present indicator (24) as employed to detect shock waves on multi-dimensional meshes. All the test cases presented in this section are governed by the Euler equations and discretized by the DG spatial discretization and the local Lax-Friedrichs type numerical flux.

2D Riemann problem
We first test the present indicator on a 2D unstructured mesh (S-1 400 ) to treat unsteady problems. The computational domain is [ −0.5, 0.5] ×[ −0.5, 0.5], S-1 400 mesh is obtained by a 400×400 structured quadrilateral mesh (mesh-size is 1 400 ), and diving each quadrilateral element into two triangles. We apply the present indicator with a fixed value (s = 0.5) and combined with the Barth-Jespersen limiter in an explicit temporal scheme to solve a 2D Riemann problem with the following initial conditions,   [20] for details of TVB-based indicators). The results obtained by different indicators are shown in Fig. 3. The results for this problem show a significant improvement with the present indicator, which leads to a better shock and contact wave resolution in the zoomed regions in Fig. 4. Furthermore, the results in Fig. 5 show the present indicator marks thinner shock and contact waves regions and fewer noises as observed in the central regions.
Readers can refer to [20] for more unsteady test cases.

2D steady flow
We now employ the present indicator on 2D unstructured meshes to treat steady problems. We apply the present indicator with a fixed fraction (θ = 0.05) and combined   results that the present indicator can detect the shock regions with high resolution and low noise. Besides, to illustrate the present indicator can lead to a stably convergent status and maintain physics properties well, residual convergence history in Fig. 8 and pressure coefficient distributions on the surface of the airfoil as compared with experiment data from [26] in Fig. 9 are provided below as well,   In addition, indication with a fixed fraction performs efficiently as combined with artificial viscosity to treat high-speed steady flows. Readers can refer to [20] for efficiency comparison and more results.

3D test cases
This part is to extend the present indicator to 3D test cases and observe its performance when applied to treat steady flows.

Inviscid flow at Ma=0.8395, α = 3.06 • around the Onera M6 wings
We use a mesh file (Tetrahedral mesh with 621508 elements) as shown in Fig. 10, DG(p2) discretization is used in spatial direction, C and β in AV term are set to 0.1 and 0.01, respectively . The convergent solution flow (density) is obtained by the present indicator (24) with a fixed fraction of 0.05, as presented in Fig. 11.(a), and the output values of the present unstructured indicator are presented in Fig. 11.(b). Residual convergence history in Fig. 12 and chordwise pressure distributions at several spanwise stations on M6 wing as compared with experiment data from [27] in Fig. 13 are provided as follows.

Inviscid flow at Ma=0.75, α = 0 • around the DLR-F4 body
We use a mesh file (Tetrahedral mesh with 550362 elements) as shown in Fig. 14, DG(p2) discretization is used in spatial direction, C and β in AV term are set to 0.01   Fig. 15.(b). Residual convergence history in Fig. 16 and chordwise pressure distributions at several spanwise stations on DLR-F4 body as compared with three groups of experiment data   Fig. 17 are provided as follows. Figure 17 shows the comparative results with experiment data. From results from CART3D and CFL3D shown in [28], we infer that the present methods by solving RANS equations or LES can obtain more perfect agreement with experimental results, which can be further developed in the future. From these two cases, one can observe the present indicator can not only detect the shock/compressing waves, but also identify the regions where characteristics compress much, such as the regions near the head and tail of the wing, which are regions that usually cause issues of convergence. As a result, employment of the present indicator in 3D can detect shock and compressing regions accurately, further leading to a stably convergent process and reduced computational costs in a steady problem.

Mesh adaptive refinement
The proposed h-adaptive method in Section 4 has been successfully applied to 2D test cases. Here, we present a simple case of transonic flow around NACA0012 airfoil to show its strength in terms of improving resolution with apparently reduced computational costs.
The initial mesh is rather coarse and is shown in Fig. 18.(a), and the initial solution (density) with low-resolution and the output of indicator are presented in Fig. 18.(b)(c). After 3 times of mesh refinement (θ=0.05), one can observe a significant improvement in terms of resolution in result (density), shown in Fig. 19.(b), and the reason can be explained from that the proposed error estimator (26) can capture shocks accurately (Fig. 19.(c)) Fig. 19 Final results of the transonic flow around the NACA0012 airfoil after 3 times of mesh refinement with θ = 0.05 and guide the computational resources (more refined mesh) focusing on the shock and compressing regions (Fig. 19.(c)).

Conclusion
We described a characteristic-featured shock wave indicator for system of conservation laws and developed it on 3D unstructured meshes. The present indicator has been successfully applied as a local error estimator in h-adaptive (mesh refinement) method to further improve resolution with less computational costs.
The numerical results demonstrated that the present indicator has excellent and robust performance on shock wave detecting in different problems. It detects shock, contact waves with low noise and high efficiency, provides an attractive alternative as an error estimator to combine with h-adaptive method to design high-resolution and high-efficiency schemes.