1 Introduction

In lots of applications, e.g., engineering analysis and design, engineers are interested in not only the whole solution but also its certain values, i.e., output functionals. For example, in computational aerodynamics, lift and drag are usually more important than the whole solution in the entire space. In such computations, engineers are interested in the accuracy of these output functionals. Considering the finite computational resources, the adaptive approach is an efficient way to improve the calculation precision of these functionals. The adaptive approach usually uses local error indicators to guide the mesh refinement. The indicators derived based on the solution features such as the large gradient in flow variables and flow discontinuities have been widely used in adaptive algorithms. Although the feature-based adaptation has had some successes, it is liable to overrefine the regions that have little effect on the output functionals^[2-3].

In order to derive reliable and efficient adaptive algorithms to produce the desired increase in the accuracy of the output functionals, adjoint-based adaptive methods have been widely devel¬oped and studied^[4-13]. Becker and Rannacher^[4-5] proposed the dual weighted residual method based on the property of Galerkin orthogonality of the finite element method. Venditti^[6] and Venditti and Darmofal^[7] derived the adaptive procedure by the local error indicators con¬structed from both the primal solution and the adjoint solution in the standard finite volume discretization. Hartmann and Honston^[8-9] and Hartmann^[10] made extensive studies on the adaptive discontinuous Galerkin (DG) method based on an adjoint approach for aerodynamic flows. Wang and Mavriplis^[11] developed the h-p adaptation approach for compressible Euler equations. Fidkowski^[12] developed anisotropic and adjoint-based mesh adaptation for the DG discretization and cut-cell meshing. However, all these studies need output-based adjoint solu¬tions as the weight in order to represent the sensitivity of the output functionals to the local primal residuals.

In shape optimization designs, the adjoint problems can be derived by two ways, i.e., the continuous approach and the discrete approach. Both of these approaches have made some suc-cesses. However, up to now, most of the adjoint-based adaptive algorithms have been developed only from the discrete adjoint approach^[3-15]. In the discrete adjoint approach, discretization is made for the primal equations firstly, and then the adjoint equations are obtained by linearizing the discrete primal problems. In this procedure, a full Jacobian matrix is needed. Therefore, if we can get it from the primal flow solver, the implementation of the adjoint equation will be much easier. However, if not, the calculation of the full Jacobian matrix will be very compli¬cated. Moreover, in the continuous adjoint approach, the adjoint partial differential equations (PDEs) are derived from the primal PDEs directly, and then these equations are discretized by the numerical schemes which may be different from the primal problem and may be inde¬pendent of the code structure of the primal problem. More importantly, we can enforce the boundary conditions more rigorously in the continuous adjoint approach. Fidkowski and Roe^[16] proposed an entropy adjoint approach, and used the continuous adjoint approach to derive the error indicators. Li et al.^[17] developed an adaptive mesh refinement algorithm, and showed the promising results based on the continuous adjoint approach in the spectral difference frame¬work. In this paper, we will provide an approach to construct the adaptive error indicator based on the continuous adjoint approach in the DG method, and show its efficiency for the error estimation.

The paper is organized as follows. In Section 2, we introduce the primal equations and the DG discretization used in this work. In Section 3, we introduce the derivation of the continuous adjoint equations and construct the adjoint-based error indicator. In Section 4, we use the continuous adjoint-based error to estimate the two-dimensional (2D) Euler equations. In Section 5, we describe the numerical implementation of our work. In Section 6, we discuss the numerical experiments to highlight the superiority of the proposed algorithm.

2 Primal equations and discretization 2.1 Governing equations

Let Ω be the computational domain in ℝ². Then, the 2D steady-state Euler equations can be written as follows:

(1)

where u=(ρ, ρu, ρv, ρE)^T is the solution vector consisting of the density ρ, V=(u, v)^T is the velocity, and E is the total energy per unit mass. The flux functional (u) is defined by

(2)

where p is the pressure. The system is closed by

(3)

where γ is the heat capacity ratio, which is 1.4 for the diatomic gas.

2.2 Standard DG discretization

Denote ={k} the subdivision of the domain Ω. Then, the standard DG discretization of Eq. (1) can be given as follows: find u_h ∈ _{h, p} such that

(4)

where

In the above equations, d is the dimension of u, ^p is the polynomial with the degree p, and (·)⁺ and (·)⁻ are the interior trace and the outer trace on ∂k, respectively. (u_h⁺, u_h^—, n) and ^b(u_h⁺, u_h^- (u_h⁺), n) are the numerical flux functionals depending on the interior and ∂Ω meeting the physical boundary conditions, respectively.

3 Adjoint-based error estimation 3.1 Continuous adjoint equations

Consider the following primal problem arising from Eq. (1):

(5)

(6)

where is a boundary operator, Γ=∂Ω, u ∈ , and is the functional space of the analytical solution. Then, the output functional which we are interested can be defined by

(7)

where i_Ω(·) and i_Γ(·) are the operators which may be nonlinear. Let u_h be the numerical approximation of u in the space _{h, p} ⊂. Then, we denote the error between u and u_h by ũ, i.e.,

For the analytical solution u u ∈ $ $ and ψ ∈ , the residual

Therefore, we can add the residual to the output functional I(u) without changing its value, i.e.,

(8)

In this way,

(9)

Therefore, the error of I(u) can be given by

(10)

Then, we linearize the right-hand side of Eq. (10) to get the first-order approximations as follows:

where the prime on (·)' means the Fréchet derivative linearized at the term in the square brackets. Thus, we have

(11)

In order to make I(u) − I(u_h) be independent of ũ, in other words, to find a suitable ψ such that the error ũ has no influence on I(u) − I(u_h), we consider the following problem: ∀ũ ∈ , find ψ such that

(12)

Therefore, ψ should satisfy

(13)

Substituting

into Eq. (13), we get

(14)

Integrating the left-hand side of the above equation by parts, we get

(15)

where the superscript ^* denotes the adjoint operation. Considering the right-hand side of Eq. (14), we can rewrite the continuous adjoint equations of Eq. (5) as follows:

(16)

(17)

3.2 Error estimation

In engineering calculation, we hope to get the best possible approximation of the output functionals with the smallest computing resources. An error estimate based on the computed values is required to evaluate how accurate these output functionals are approximated. From prior experience, it is very hard to estimate the error. Using the adjoint problem, we can obtain the general result.

Theorem 3.1 Let u be the analytical solution and u_h be the numerical solution of the primal equation (5). Let ψ denote the analytical solution of the adjoint equation (16). Then, the error between I (u) and I (u_h) can be estimated by

(18)

where

(19)

From this theorem, we can conclude that the primal residual can be related to the output functional error via the adjoint variable as a weight. Then, we take the absolute value of ηk to construct the local error indicator for adaptive algorithms, namely,

(20)

4 Application on Euler equations

4.1 Some basic notations and output functionals

In this section, firstly, we define the flux Jacobian matrix as follows:

Secondly, we define the boundary flux as follows:

where n=(n₁, n₂) is the normal vector of the boundary, and u_n is the normal velocity, i.e., u_n=un₁ + vn₂. The pressure induced force coefficients, such as the drag coefficient, are the most important output functionals in inviscid compressible flows. Here, the drag coefficient is given by

(21)

where Γ=∂Ω, and Γ_W is the solid wall of Γ. where l is the reference length. ϕ=(cos(θ), sin(θ))^T, where θ denotes the angle of attack. The subscript ∞ means the far-field quantity.

4.2 Adjoint equation and boundary conditions

In this section, based on the drag coefficient, we introduce how to derive the continuous adjoint Euler equations and its boundary conditions. The approaches for other output functionals are similar. From Eq. (21), we note that the drag coefficient does not depend on the domain terms, i.e., i_Ω(u)=0 in Eq. (7). Therefore, the adjoint equation (16) can be given by

(22)

where ψ=(ψ₁, ψ₂, ψ₃, ψ₄) is the adjoint solution vector.

In inviscid compressible flows, on the solid wall Γ_W, the boundary condition imposes u_n=0. Therefore, n. (u)=p(0，n₁，n₂, 0)^T on Γ_W. Recalling the general functional of the boundary condition (17), the adjoint boundary condition on Γ_W under this output functional can be given by

(23)

i.e.,

(24)

On a far-field boundary ∂Ω_∞, recall the definition of the drag coefficient (21). We notice that the value of this output functional does not depend on the far-field. Therefore, simply, we set

(25)

Finally, we obtain the adjoint problem as follows:

(26)

5 Numerical implementation 5.1 Solving primal equations

Firstly, we add a pseudo time term to Eq.(1), i.e., we rewrite the steady-state Euler equations as follows:

(27)

The spatial discretization by the DG method leads to a system of ordinary differential equations as follows:

(28)

where is the mass matrix, and R is the right-hand residual. In this work, we are interested in stationary problems. Therefore, the backward Euler scheme is used in the temporal discretization. After temporal discretization and linearizing Eq. (28) in time, we can get a system of linear equations as follows:

(29)

where Δt is the time increment, and

Then, the linear system can be solved by the generalized minimal residual (GMRES) method with ILU0 preconditioning.

5.2 Solving adjoint equations

In Eq.(19), the error indicator involves the analytical adjoint solution ψ which is unknown. Therefore, we need to make an approximation to it. Let ψ_h denote the approximation. We note that ψ_h must not be approximated by the same finite element space V_{h, p}. Otherwise, the error indicator will be identically equal to zero due to the Galerkin orthogonality. There are several approaches to compute the numerical approximation ψ_h^[14-15]. In this work, the adjoint equations are approximated on the same mesh employed for the primal equations, but the order of the finite element space is increased from p to p + 1.

The numerical implementation of the adjoint equations is a little different from the primal equations. From Eq. (26), we notice that (A₁^T, A₂^T) is independent of the adjoint variable ψ. Therefore, the adjoint equations are linear. Therefore, we can solve this problem by iteration or the direct linear system solver directly after spatial discretization.

We should remark that there are two main differences between the DG discretization of the primal problem and the DG discretization of the adjoin problem. The first one is the difference between the numerical fluxes. Here, we write Eq. (1) in the following quasi-linear form:

i.e., (A₁, A₂) · ∇u=0 in Ω. We note that Eq. (22) and the above equation have the same eigenvalues. Therefore, the local Lax-Friedrichs numerical flux for the adjoint equation can be given by

(30)

where α is the largest eigenvalue which can be calculated from the primal problem directly. The second difference lies in the boundary conditions. As we have discussed in Section 4 on the far-field boundary ∂Ω_∞, we simply set ψ_h⁻=0. On the solid wall boundary, there are a set of boundary values of ψ_h satisfying the boundary conditions in Eq. (26). Here, we enforce the solid wall boundary condition as follows:

(31)

(32)

(33)

(34)

which guarantees the average value

satisfying the adjoint solid wall boundary condition, i.e., ψ₂n₁ + ψ₃n₂=

5.3 Adjoint-based adaptive mesh refinement algorithm

Utilize the local error indicator (20) by an adaptive mesh refinement algorithm shown in Fig. 1. In each refinement cycle, the cells with a large local error indicator η_k^ind will be refined, and the cells with small η_k^ind ind will be coarsened. Following the process, we will get a mesh on which the local error indicator η_k^ind is equidistribution. For more details of the mesh adaptation strategies, we refer readers to Ref.[8].

6 Numerical results

We show some numerical tests to discuss the advantages of our adaptive mesh refinement algorithm, comparing with the traditional algorithms based on the feature-based or residual-based error indicators. Here, firstly, let us introduce the feature-based error indicator η_f^eature=lg(1 + | ∇ρk |) and the residual-based error indicator^[18]

where ∇ρk is the gradient of the density, and

In the following numerical tests, the primal problem will be approximated in _{h, 1}, and the adjoint problem will be approximated in _{h, 2}.

6.1 NACA-0012 airfoil

In the first numerical experiment, we consider an NACA-0012 airfoil in an inviscid compressible flow with Ma_∞=0.4 and θ=0. The output functional considered here is the drag coefficient. Therefore, the theoretical value of the output functional is zero.

Figure 2 shows the initial mesh of this test. The meshes obtained by use of both the adjoint error indicator and the feature-based error indicator are shown in Fig. 3. From these figures, we can see that most of the refinement is enforced near the leading and trailing edges of the airfoil. Furthermore, we find that the mesh obtained by our adjoint-based error indicator is also refined around the surfaces of the airfoil.

Figure 4 shows the output functional error convergence under different refinement approaches. From these results, we can see that the error of the output functional estimated on the mesh obtained by use of the adjoint-based error indicator is the smallest for a given number of cells. In addition, we can observe that, at a given error, the number of the mesh cells needed by the adjoint-based adaptive approach is an order of magnitude less than the uniform refinement approach.

Figure 5 shows the distribution of the adjoint-based error indicator during the refinement cycle. We find that the errors are larger around the surfaces of the airfoil on the initial mesh.

The cells with larger errors are refined, which explains why the cells near the airfoil surfaces are refined. Also, we notice that with the repetition of the mesh refinement, the errors in all cells are decreasing. Consequently, using our adaptive algorithm, we can get a mesh, on which the error is equally distributed.

6.2 Smooth bump

In the second numerical experiment, we consider the inviscid flow over a smooth bump. The computational domain is given by Ω={(x, y) ∈ [−1.5, 1.5] × [0, 0.8], y > 0.0625 e^−25x2}. At the inlet, Ma_∞=0.5, and θ=0. There should be no entropy production because the flow is smooth enough. Therefore, we construct the following output functional:

where p is the pressure, ρ is the density, and the subscript ∞ means the freestream value. At the left-and right-sides, we apply the subsonic inflow boundary condition and the subsonic outflow boundary condition, respectively. At the top and the bottom, the slip wall conditions are applied. The theoretical value of this output functional is zero.

Figure 6 shows the initial mesh for this problem, and Fig. 7 shows the pressure contour of the primal Euler equations. We can see that the solution is axisymmetric. From Fig. 8, we can see that most of the refinements are near the bump. Both of the two adaptive approaches construct almost the same meshes, except that the refinement is denser around the bump in the adjoint-based approach.

In Fig. 9, we show the convergence of the output functional error. Some observation can be made. Firstly, both the adjoint-based and residual-based adaptive approaches show better convergence than the uniform refinement. Indeed, after several refinements, the error computed on the subsequent adaptive mesh is at least two orders of magnitude smaller comparing with the error computed on the uniform refinement mesh. Secondly, the adjoint-based refinement gives a better accuracy than the residual-based refinement for a given number of cells.

7 Conclusions

We have introduced an adjoint-based adaptive mesh refinement algorithm for inviscid compressible flows. Firstly, we show the derivation of a continuous adjoint problem for the Euler equations. Secondly, we introduce the error estimation for the given output functionals. Thirdly, we discuss the numerical implementation of the DG method for an adjoint problem. Then, we develop an adaptive mesh refinement algorithm based on the continuous adjoint approach. The numerical tests show that the meshes constructed by this adaptive algorithm are much more economical for computing the values of the output functionals than the meshes constructed by the traditional error indicators. We conclude that the adaptive approach presented in this paper is very promising in improving the accuracy of the output functionals within given computational resources.