Smoothed truncation error in functional error estimation and correction using adjoint methods in an unstructured finite volume solver

Computers and Fluids 140 (2016) 406–421

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

Smoothed truncation error in functional error estimation and correction using adjoint methods in an unstructured finite volume solver Mahkame Sharbatdar1,∗, Alireza Jalali2, Carl Ollivier-Gooch3 Department of Mechanical Engineering, The University of British Columbia, Vancouver, BC, V6T 1Z4, Canada

a r t i c l e

i n f o

Article history: Received 22 March 2016 Revised 9 September 2016 Accepted 19 October 2016 Available online 19 October 2016 Keywords: Output functional correction Adjoint methods Smooth estimation of truncation error Unstructured meshes Finite volume methods

a b s t r a c t The goal of this paper is to develop a reliable truncation error estimator for finite-volume schemes and to use this truncation error estimate to correct the output functional of interest. The rough modes in the truncation error for unstructured mesh are dominant and if p-truncation error estimation is used, functional correction has a poor performance. So, we are trying to obtain a smooth estimate of the truncation error to improve the performance of output error estimation. The correction term is based on the truncation error and the adjoint solution. Both discrete and continuous adjoint solutions are used for correcting the functional. Our results for a variety of governing equations suggest that the smoothing scheme can improve the correction process significantly. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction In the field of computational aerodynamics and fluid dynamics, lift and drag are output functionals of interest and the desire for efficient computational methods producing reliable and accurate lift and drag values drives algorithm research in the field. The adjoint method has played an important role in this context because of the great flexibility it offers with regard to the physics model and to the definition of output functionals. The history of the use of adjoint equations in fluid dynamics design goes back to the work by Pironneau [35] and particularly in the field of computational aerodynamics design to the work by Jameson [16]. Since then, adjoint methods have been used for design applications for both internal and external flows [15,16,18,19,23,36]. The adjoint theory was first presented in the context of linear algebra by using the algebraic equations obtained from the discretization of the original problem. This is the basis for the discrete adjoint approach. The continuous adjoint approach, on the other hand, is formulated based on the adjoint PDE which is discretized and solved independently [11].

∗

Corresponding author. E-mail addresses: [email protected] (M. Sharbatdar), interchange.ubc.ca (A. Jalali), [email protected] (C. Ollivier-Gooch). 1 PhD Candidate. 2 PhD Candidate. 3 Professor. http://dx.doi.org/10.1016/j.compfluid.2016.10.019 0045-7930/© 2016 Elsevier Ltd. All rights reserved.

arjalali@

The adjoint problem also plays a key role in estimating and reducing the error in output functionals. Within the context of finite element methods, output functional correction has been outlined by Becker and Rannacher [5] and Larson and Barth [20] based on structural finite element methods [2]. The adjoint-based error correction technique developed by Pierce and Giles [33] extends the inherent super-convergence properties of finite element methods to cover numerical results obtained from finite difference and finite volume methods without natural super-convergence properties. Moreover, the technique can be used to improve the accuracy of super-convergent functionals obtained from finite element methods by constructing smoother, higher-order interpolants of the primal and adjoint solutions [33]. This method for output error estimation and correction for two-dimensional inviscid flows was applied by Venditti and Darmofal [41] for a second-order finite volume discretization. For numerical efficiency, computing the functional value to a higher order of accuracy than the primal solution is advantageous. A super-convergent estimate of the functional may be obtained by computing the leading error term in the original functional estimate and using this as a correction. Pierce and Giles [34] showed that the error in the functional value based on the reconstructed primal solution can be expressed as a function of the truncation error, implying that truncation error estimation is required for output error estimation. Truncation error is defined as the difference between the continuous PDE and the finite discretized equation. Historically,

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

truncation error analysis for structured mesh schemes is a routine application of Taylor series analysis, while analysis of unstructured mesh schemes has lagged behind. Recent work on analysis of unstructured mesh schemes includes accuracy of discretization schemes on irregular grids [8], error comparisons for cell-centered and vertex-centered discretizations [9], and Taylor-based accuracy assessment of cell-centered discretizations [14]. Ollivier-Gooch and Van-Altena [31] performed Taylor series truncation error analysis for the Laplacian on regular triangular meshes, and this was extended to general stencils [13,14]. The behavior of truncation error on unstructured meshes is completely different from structured meshes. Several researchers, including Diskin and Thomas [7,9], and Jalali and Ollivier-Gooch [14], have demonstrated that the truncation error for unstructured finite volume schemes is asymptotically larger than the discretization error, which in turn is typically of the same order as the solution approximation error in the scheme. This behavior is in contrast with the structured mesh case for which the truncation error has the same asymptotic order of accuracy as the discretization error. Another feature of the truncation error for unstructured mesh schemes is its noisy appearance, caused by the discontinuous jump of the coefficients of the terms in the Taylor series expansion of the error from one control volume to another; this is in contrast with structured mesh schemes where the truncation error is smooth. These two features of the truncation error for unstructured mesh finite volume schemes, non-smoothness and large magnitude, are related to each other by the eigensystem of the discrete problem as shown by Ollivier-Gooch and Roy [30]. Sharbatdar and Ollivier-Gooch [39] have shown by eigenanalysis of the truncation error that the rough modes dominate the unstructured mesh truncation error. The dominance of the rough modes causes difficulties in output error estimation. We use a higher-order flux integral as an estimation of the truncation error, the p−TE method, and the rough modes are responsible for the difficulty in accurately estimating discrete truncation error. We try to develop a smooth estimation of the truncation error which can be exploited in output functional correction. Venditti and Darmofal presented an error estimation strategy based on the adjoint formulation for estimating errors in functional outputs for one-dimensional problems and their error estimation procedure was applied to a standard, second-order finite volume discretization [40]. The goal of this paper is to investigate continuous truncation error estimates by producing a continuous approximation to the finite volume discrete solution and estimating truncation error as the higher-order residual of the PDE for this continuous function. The truncation error can be used to find the leading error term in the functional and we compare the performance of continuous and discrete adjoint approaches in the correction procedure. We begin by briefly describing our higher order finite volume flow solver followed by adjoint problem definition and implementation in Section 3. The output functional calculation process and error correction calculation is described in Section 4. We will show that the correction term is a function of the truncation error and our method for estimating and improving the truncation error is illustrated in Section 5. Several test cases with different governing equations and physical behaviors are shown in Section 6. We show the capability of the correction scheme we developed for unstructured mesh finite volume scheme for both scalar equations and system of equations; advection and Poisson as examples of scalar equations and Euler and Navier–Stokes for system of equations.

2. Higher order finite volume flow solver To discretize the flow equations using the finite volume method, the governing equations should be recast in fully conser-

407

Fig. 1. Schematic illustration of Gauss quadrature points for third and fourth-order .

vative form as:

∂ U + ∇ · F = 0 ∂t

(1)

where U denotes the solution vector and F is the flux vector. Integrating Eq. (1) over an arbitrary control volume and using the divergence theorem gives the 2D finite volume formulation of the governing equations in the form of a volume and a surface integral

¨ CV

∂ U dA + ∂t

‹ CS

F · nˆ ds = 0

(2)

where nˆ represents the outward unit normal vector, CV and CS are the control volume and control surface (the boundary of the control volume), respectively, and ds is the infinitesimal. Assuming that the discretized physical domain does not change in time, the time derivative can be brought out from the integral in Eq. (2), giving an evolution equation for the vector Ui , the control volume average solution,

1 dU i =− dt ACVi

˛

CSi

F · nˆ ds

(3)

The right hand side is called the flux integral or residual and represents the spatial discretization of the same control volume. To compute the flux integral for each control volume, the numerical flux is computed at Gauss quadrature points based on the reconstructed solution. For second-order, one quadrature point is required and for third and fourth order flux integration, two quadrature points per edge are necessary. Fig. 1 shows the higher-order quadrature points for the flux integration schematically. To obtain a single flux from two different values at each Gauss point, we use Roe’s scheme [38] for inviscid fluxes; for viscous fluxes, we average the gradients and add a jump term [27]. The control volume flux integral in Eq. (3) is approximated as the summation of flux integrals over the edges and can be re-written as:

  dU i = −Ri U i dt

(4)

 

where Ri U i is called the residual.   For steady-state problems, R U¯ i = 0. One could solve the nonlinear system of equations by the direct application of Newton’s methods for steady-state problems. However, Newton’s method will diverge if the initial guess is too far from the real solution. As a result, the non-linear system is augmented by a damping term which mimics the time derivative in the original time-dependent equations and prevents the evolution of non-physical solution at each iteration. This is called the implicit pseudo time-stepping method [22]. For implicit Euler integration of the discretized equation in time, the next time level solution for both the space and the time discretizations are used. Assuming that the solution avern age vector at the current time level n is denoted by U i , both sides of Eq. (4) should be evaluated at the next time level n + 1. After linearization, we get



I

t

+

 ∂R n δUi = −Ri (U i ), ∂U

n+1

Ui

n

= U i + δUi

(5)

408

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

where ∂∂UR is the Jacobian matrix resulting from flux integral linearization. Eq. (5) is a large linear system of equations which needs to be solved at each time step to obtain an update for the vector of unknowns. The global Jacobian matrix can be represented explicitly as:

∂ R ∂ F luxInt ∂ F lux ∂ RecSol ∂ RecCoe f = ∂ F lux ∂ RecSol ∂ RecCoe f ∂U ∂U

(6)

The last term is found by using the pseudoinverse of the reconRecSol struction matrix, ∂∂RecCoe is a geometric term dependent on the f

´ where (U, V )D = D U T V dv is an integral over the domain and ´ (U, V )∂ D = ∂ D U T V dA is an integral over the boundary of the domain. The objective in the analytic approach is to convert the primal problem using integration by parts into an equivalent adjoint problem

Determine J = (Z, f )D + (C ∗ Z, e )∂ D

given that L∗ Z = g ∗

B Z=h

in D on ∂ D

L∗

(11) B∗

F lux location of the Gauss point, ∂∂RecSol is based on the discretization ∂ F luxInt scheme and ∂ F lux is computed using the appropriate Gauss integration weight and the length of edges. More details concerning the higher order finite volume flow solver that we are using are given in Refs. [22,24,29].

where is the linear PDE which is adjoint to L. B and are boundary condition operators for the primal and adjoint problems, respectively, and C and C∗ are possible differential weight boundary operators used in the functionals; these four operators may have different dimensions on different parts of the boundary. The two forms of the problem are equivalent provided that

3. Continuous and discrete adjoint methods

J = (Z, LU )D + (C ∗ Z, BU )∂ D = (L∗ Z, U )D + (B∗ Z, CU )∂ D

The adjoint theory in the context of linear algebra is the basis for the discrete adjoint approach. In this approach, the algebraic equations arising from discretization of the original fluid dynamic equations are used. Alternatively, the adjoint PDE can be formulated and discretized independently, which is the basis of the continuous adjoint formulation. These two different adjoint approaches are described briefly in this section. 3.1. Discrete adjoint Consider a partial differential equation arising from a finite volume discretization of the original fluid dynamic equations which is discretized and written as an algebraic equation:

 

Rh U h = 0



(7)



(12)

which is the adjoint consistency condition. The left hand side is the output functional based on the adjoint solution and the right hand side is the output functional based on the primal problem. The continuous adjoint has been used for optimizing wing shape [17] and correcting the output functional [3]. Comparing the discrete adjoint to the continuous adjoint reveals that although the discrete adjoint solves a system of equations based on the transpose of the Jacobian matrix and the derivative of the functional with respect to the solution, the continuous adjoint is an independent PDE which needs to be discretized and solved by itself. It implies that boundary conditions should be set for the continuous adjoint problem and those boundary conditions are dependent on the functional. In Section 6, boundary conditions for different problems based on corresponding functionals are discussed.

For a scalar output, Jh U h such as lift or drag, the associated discrete adjoint vector, Zh , must satisfy the discrete adjoint equation:

4. Output functional correction



Pierce and Giles have developed a method to estimate error in the functional [34] and an overview of the method is discussed here. Considering Eq. (12) and for given approximate reconstructed solutions Uh and Zh :

∂ Rh ∂U h

T



∂ Jh Zh + ∂U h

T =0

(8)

A detailed description of the discrete adjoint formulation is given by Fidkowski and Darmofal [10]. As discussed in Section 2, the implicit Jacobian matrix can be computed in our solver and hence the ∂R transpose of the Jacobian matrix, ( h )T , can be calculated easily. ∂U h

The discrete adjoint solution can be used for mesh adaptation and optimization [3,10,25,26]. If we assume that both the equations and the functional J have been linearized, the discrete approach can be described as a mapping from the original problem J = (U, g) given that AU = f into an equivalent adjoint problem J = (Z, f ) given that AT Z = g. The inner product is a vector dot product (V , U ) = V T U and the equivalence of the two problems is easily proved to be





(Z, f ) = (Z, AU ) ≡ AT Z, U = (g, U )

(9)

Note that the inhomogeneous term f in the discrete equations in the primal problem, enters the functional in the adjoint problem, and correspondingly the inhomogeneous term g in the adjoint problem comes from the functional of the primal problem. 3.2. Continuous adjoint

Determine J = (U, g)D + (CU, h )∂ D BU = e

in D in ∂ D

= (L∗ Zh , Uh )D + (L∗ Zh , U − Uh )D + (g − L∗ Zh , U − Uh )D +(g − L∗ Zh , Uh )D

(B∗ Zh , CUh )∂ D + (B∗ Zh , C (U − Uh ))∂ D + (h − B∗ Zh , C (U − Uh ))∂ D +(h − B∗ Zh , CUh )∂ D (13) By using integration by parts, we can re-write the above equation as

J = (L∗ Zh , Uh )D + (B∗ Zh , CUh )∂ D

+(Zh , L(U − Uh ) )D + (C ∗ Zh , B(U − Uh ) )∂ D +

(Z − Zh , L(U − Uh ))D + (C ∗ (Z − Zh ), B(U − Uh ))∂ D

(10)

(14)

The two terms in the first line are the influence of bulk and boundary integrals in the discrete functional, the terms in second line represent computable adjoint error estimates and the last two terms in the third line are the higher order remaining errors which are negligible compared to other terms. Accordingly, the estimated error in calculating the functional is:

δ J = (Zh , f − LUh )D + (C ∗ Zh , e − BUh )∂ D = (Zh , τ )D + (C ∗ Zh , e − BUh )∂ D

The continuous primal problem can be posed as

given that LU = f

J = (g, Uh + (U − Uh ) )D + (h, C (Uh + (U − Uh ) ) )∂ D

(15)

where τ is truncation error. Hence, an estimate of the truncation error is needed to find the correction term for the output functional.

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

409

5. Truncation error estimation and smoothness

5.1. Discrete p-Truncation error estimation One of the methods that is generally applicable for estimating the truncation error on unstructured meshes is the p-TE method, which uses higher-order flux integrals. We can solve the original discrete problem to order p, where the residual computed to this accuracy is necessarily zero:





Lh,p Uh,p = 0

(16)

Lh, p is a discrete approximation with order of accuracy p to a partial differential operator on a mesh with characteristic size of h. For estimating the truncation error, a higher order discretization is used to find the leading order terms in the truncation error estimate:



Lh,p+1 (Uh,p ) = τh;p→ p+1 Uh,p



(17)

Consider that the p-order operator is applied to the higher-order solution:







Lh;p Uh;p+1 = Lh;p Uh;p + δU



(18)

where δ U is the difference between the two solutions, δU = Uh,p+1 − Uh,p . Eq. (18) can be expanded by using the linearity of the operator as:









Lh;p Uh;p+1 = Lh;p+1 − Lh;p+1 − Lh;p







Uh;p + δU







= Lh;p+1 Uh;p+1 − Lh;p+1 − Lh;p Uh;p + δU





















= −Lh;p+1 Uh;p + Lh;p Uh;p − Lh;p+1 (δU ) + Lh;p (δU ) = −Lh;p+1 Uh;p + 0 − Lh;p+1 (δU ) + Lh;p (δU )

(19)

For a structured mesh, as a result of symmetry and smoothness in the mesh, the problem is well-behaved and the last two terms in the above equation are of order of p + 1 and p, respectively. Hence, for a structured mesh,





















Lh;p Uh;p+1 = −Lh;p+1 Uh;p − O h p+1 + O(h p ) ⇒ Lh;p Uh;p+1 = −Lh;p+1 Uh;p + H.O.T

-4

10

(20)

where the left hand side is the leading-order term in the exact truncation in the sense of being the truncation error associated with the leading-order term in the discretization error, and the right-hand side is a higher-order flux integral based on the pthorder solution. Estimating the truncation error by Eq. (20) is called

2.00

10

-6

4.70

J str Corr Jstr J unstr Corr Junstr

-8

40

60

80

(#CV)1/2

100

120 140

Fig. 2. Comparison of output error estimation for structured and unstructured meshes .

the p−TE method. However, truncation error for the unstructured mesh finite volume method is asymptotically larger than the discretization error and contains a large amount of noise such that there is no smoothness in the truncation error. The rough modes in an eigenanalysis based on the Jacobian matrix are the dominant modes in the truncation error and numerical experiments show that by applying the linear operator on δ U, we lose orders in computing the last two terms in Eq. (19). Details on this eigenanalysis and the dominance of the rough modes of the truncation error on unstructured meshes is described by Sharbatdar and Ollivier-Gooch [39]. The effect of the different behavior of truncation error estimation for structured and unstructured meshes on output error estimation can be shown for a Poisson test case on a (1 × 1) domain. A Poisson problem with an exact solution is considered:

∂ 2u ∂ 2u π3 + = − sin (π x ) sin (π y ) 4 ∂ x2 ∂ y2 π ⇒ uexact = sin (π x ) sin (π y ) 8

and the boundary conditions are set as zero to match the exact solution. The weight function which is also the source term of the adjoint problem, is set so that the exact output functional is one:

π5

x(1 − x )y(1 − y ) 2 ˆ 1ˆ 1 π J = (u, g) = sin (π x ) sin (π y ) 0 0 8 g=



= 0 − Lh;p+1 − Lh;p Uh;p − Lh;p+1 − Lh;p (δU )



10

Error

There are several options for estimating the truncation error. If the exact solution is available, the exact continuous solution can be used in the discrete residual which provides us the exact truncation error. This is useful for comparing other estimates with, but is not applicable to general complex problems. Instead of using the exact solution, the discretization on a finer mesh can be used, this method is called h-TE method. Darmofal et al. [41], Nemec and Aftomis [25] and Hartman [21] used this method for finite element adjoint-based adaptation. The requirement of the residual on the fine mesh implies more computational costs in terms of time and memory for h-TE method. The method that we are using in this paper is called the p-TE method, in which a higher order discretization is used for estimating the truncation error. This approach has been used for truncation error estimation for structured meshes [33], but behaves in a completely different way on unstructured meshes. In this section, the p−TE method is described and we will show how different and inaccurate it is for unstructured meshes. The approach we developed for improving this estimate is then described.

×

π5 2

x(1 − x )y(1 − y )dxdy = 1

Fig. 2 compares the original functional based on the second order solution, black lines, and the corrected functionals using fourth order flux integrals to estimate the truncation error, red lines, for both structured and unstructured meshes of approximately the same number of control volumes. Although using a fourth-order correction term on unstructured meshes reduces the error in calculating the functional, the convergence rate does not improve with mesh refinement. However, using the structured mesh, the rate of convergence is much faster and even better than quartic convergence. This totally different behavior for structured and unstructured mesh is the result of different smoothness properties of the truncation error on these two sets of meshes.

410

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

Fig. 3. Argyris reference element.

5.2. Estimate of the truncation error from continuous data Because the disappointing results shown in Fig. 2 are associated with lack of smoothness of the truncation error, we seek to develop a smooth truncation error estimate that can be computed using quantities already available in a typical finite volume flow solver. Specifically, we have chosen to compute a C1 interpolation that approximates the reconstructed finite volume solution. We then use this interpolation to evaluate the analytic fluxes and integrate them to compute a truncation error estimate. To compute the C1 interpolation, we use the Argyris element because of its high rate of spatial convergence (O(h6 ) in the L2 norm for sufficiently smooth problems) and simplicity of implementation compared to other elements. The Argyris triangle [4,6] is based on the fifth order space of quintic polynomials over a triangle. This element produces C2 continuity at the vertices and full C1 continuity along the edges between elements of a triangulation which suffices for our case. Quintic polynomials in R2 are a 21-dimensional space, and can be evaluated by using six pieces of data per vertex and one per edge. The vertex degrees of freedom are the solution, two first derivatives to specify the gradient, and three second derivatives to specify the unique components of the (symmetric) Hessian matrix; the edge degree of freedom is the normal derivative as shown in Fig. 3. To compute these values, we first reconstruct the finite volume solution to higher (third or fourth) order. For each edge or vertex, the reconstruction gives multiple values for the solution and its derivatives, one per cell. For a pth order reconstruction, these solution values will differ only by O(hp ) for smooth data, with derivatives losing one order. We take an arithmetic average of these values as our input to the interpolation. The Argyris element has 21 basis functions, each of which is non-zero for exactly one of the 21 degrees of freedom. The solution is then represented as

   ∂ u  ∂ u  ∂ 2 u  u = u0 ϕ0 + ϕ + ϕ + ϕ ∂ξ 0 1 ∂η 0 2 ∂ξ 2 0 3   ∂ 2 u  ∂ 2 u  + ϕ + ϕ ∂ ξ ∂ η 0 4 ∂η2 0 5    ∂ u  ∂ u  ∂ 2 u  +u1 ϕ6 + ϕ + ϕ + ϕ ∂ξ 1 7 ∂η 1 8 ∂ξ 2 1 9   ∂ 2 u  ∂ 2 u  + ϕ + ϕ 10 ∂ ξ ∂ η 1 ∂η2 1 11    ∂ u  ∂ u  ∂ 2 u  +u2 ϕ12 + ϕ + ϕ + ϕ ∂ξ 2 13 ∂η 2 14 ∂ξ 2 2 15

Fig. 4. Comparison of L2 -norm of the exact truncation error and the estimation of the truncation error calculated by 4th-order flux integral based on smoothed and non-smoothed solution.

  ∂ 2 u  ∂ 2 u  + ϕ + ϕ (21) ∂ ξ ∂ η 2 16 ∂η2 2 17       ∂u ∂u ∂u ∂u + n + n  ϕ + n + n  ϕ ∂ξ ξ ∂η η  1 ,0 18 ∂ξ ξ ∂η η  1 , 1 19 (2 ) (2 2)    ∂u ∂u + nξ + nη  ϕ20 ∂ξ ∂η (0, 12 )

where

ϕi (x, y ) = c + cx x + cy y + cx2 x2 + cxy xy + cy2 y2 + cx3 x3 + cx2 y x2 y + cxy2 xy2 + cy3 y3 + cx4 x4 + cx3 y x3 y + cx2 y2 x2 y2 + cxy3 xy3 + cy4 y4 + cx5 x5 + cx4 y x4 y + cx3 y2 x3 y2 + cx2 y3 x2 y3 + cxy4 xy4 + cy5 y5 (22) More details about calculating the 21 coefficients for each of these basis functions can be found in Appendix A. By using this C1 interpolated solution, a single value for the solution and gradients is obtained at the quadrature points of two neighboring cells, Figure 1. We compute a p−TE estimate which is the higher order flux integral calculated by using this smooth flux vector at each quadrature point. Although the typical p−TE estimate uses higher-order flux integral based on the discrete lower-order solution, the new developed estimate of the truncation error uses the higher-order flux integral based on C1 interpolation of the lower-order discrete solution. Fig. 4 shows the L2 −norm of this estimation of the truncation error calculated by fourth-order flux integral based on the secondorder solution for the Poisson problem. This estimation of the truncation error is much closer to the exact truncation error and the convergence rate is much faster than the estimation of the truncation error based on the non-smoothed solution. In calculating the flux integral for the Poisson problem, one order of accuracy is lost in calculating the flux (the gradients) and another order is lost in calculating the flux integral. As a result, the fourth-order flux integration based on the smoothed second-order solution is secondorder. On the other hand, the estimation of the truncation error

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

based on the non-smoothed solution is not converging with the expected rate. The estimate of the truncation error based on the smoothed solution is multiplied by the adjoint solution to find the correction term added to the output functional to have a more accurate functional which converges to the exact value with the higher-order rate. It should be noted that the correction term of Eq. (15) is the summation over all control volumes of the truncation error calculated by higher-order flux integral based on lower-order solution multiplied by the adjoint solution of the same order as the primal problem

δ J = (Zh , τ )D =

 ˛

Fpp+1



· nˆ ds Z¯ p

CV

The higher-order flux integration can be replaced by p + 2 as well. 6. Results and discussion We conducted several experiments for different governing equations with different physical behavior, for both convective and diffusive fluxes, to show the effectiveness of our smoothing scheme in correcting the output functional, compared to functionals corrected by the non-smooth p−TE in unstructured meshes. Both scalar and system problems are considered. For each case, the governing equation and the domain of the problem are introduced and then correction based on non-smooth and smooth truncation error are compared for both discrete and continuous adjoint problems. To verify our approach to improving the functional, we begin with problems with exact solutions so that we can calculate the exact functional to compare the magnitude of the error and also the convergence rate of the corrected functional with mesh refinement. Since this technique is mathematical in origin rather than physical, applying it to other problems, whether they have exact solutions or not, follows exactly the same procedure. We are focusing on problems with exact solutions only to have an exact reference value for comparison purposes. Furthermore, since adjoint consistency is a key element of discrete adjoint methods, adjoint consistency should be verified. Instead of showing adjoint consistency for each problem separately, we will show it only for the Navier–Stokes equations which is a real flow with both convective and diffusive fluxes. 6.1. Scalar equations We start by examining convective and diffusive fluxes separately using the advection and Poisson problems, respectively. Hence, two model problems are defined by the method of manufactured solutions [37]. 6.1.1. Advection equation Consider the 2D steady linear advection equation with zero source term:

∇ · (bu ) = 0 in D u = uin on ∂ Di

(23)

where b is the velocity vector and ui is the solution at the inlet ∂ Di . The left hand side of Eq. (23) is multiplied by z, integrated over the domain D and integrated by parts to obtain the continuous adjoint equation:

    (∇ · (bu ), z )D + u, −b · nˆ z ∂ D = (u, −b · ∇ z )D + u, b · nˆ z ∂ Do i

(24) where nˆ is the outward normal vector. Comparing to Eq. (12), we end up with Lu = ∇ · (bu ) in D and

Bu = u,

Cu = 0

on ∂ Di

Bu = 0,

Cu = u

411

on ∂ Do

for the primal problem and L∗ z = −b · ∇ z in D and

B∗ z = 0, ∗

on ∂ Di

C ∗ z = −b · nˆ z

B z = b · nˆ z,

on ∂ Do

∗

C z=0

for the continuous adjoint problem. The output functional is defined as

ˆ

J=

D

ˆ =

D

(∇ · (bu ))zdA + (−b · ∇ z )udA +

ˆ

ˆ

∂ Di ∂ Do





uin −b · nˆ z ds





b · nˆ z uds

(25)

Setting the velocity vector as b = (1, 0 ), the adjoint problem would be an advection problem with reverse velocity of the primal problem. The inlet velocity of the primal problem is set as sin (π y), so the primal problem is

∂u = 0 ⇒ uexact = sin (π y ) ∂x The problem domain considered is a (5 × 1) rectangular channel. The output functional is defined as

J = (u, g)D + (Cu, h )∂ D0 = 0 + (u, h )∂ D0 ˆ 1 = sin (π y ) sin (π y )dy = 0.5 0

implying that the adjoint problem boundary condition should be sin (π y). The functional is defined intentionally such that the adjoint problem has exactly the same distribution as the primal problem. Both continuous and discrete adjoint solutions are depicted in Fig. 5. Although both solutions are qualitatively similar to the desired solution, the continuous adjoint is smoother and more accurate, as expected, as a consequence of solving the continuous adjoint PDE instead of solving for the transpose of the Jacobian arising from the discrete solution. This increased accuracy is reflected in the correction process as well. Before showing the numerical results, the nomenclature we are using should be described; J is the functional, CJ is the corrected functional and C’s superscript represents whether the correction is based on truncation error using the non-smoothed solution, Cn , or is based on the truncation error using the smoothed solution, Cs . J’s subscript shows the primal solution order of accuracy and the superscript shows the order of reconstruction and flux integration used as an estimation of the truncation error. The error in the corrected functional is depicted in Fig. 6 using both continuous and discrete adjoint solutions. Comparing the functional based on second-order solution with the corrected functional based on the non-smoothed solution shows poor performance of the correction in terms of error magnitude and also asymptotic convergence rate for both continuous and discrete adjoint solutions. On the other hand, correcting the functional by using the truncation error based on the smoothed solution, C s J23 , has smaller error compared to the second-order functional and the rate of convergence is 4.06 and 3.31 for continuous and discrete adjoint, respectively. In this case, we have super-convergence with the continuous adjoint and the convergence is quartic instead of cubic. Using fourth order flux integration for estimating the truncation error based on the smoothed second-order solution, C s J24 , further reduces the error. Comparing the error magnitudes of these corrected solutions based on second-order solution with the error of the third-order solution, reveals that this correction scheme based on the smoothed solution is even more accurate than the higherorder functional. Instead of using the second-order functional, the third-order functional may be corrected based on fourth-order flux integration, C s J34 . Similarly, the convergence is quartic for the corrected solution although the original functional is cubic.

412

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

Fig. 5. Adjoint solutions for advection problem.

Fig. 6. Convergence history for the advection problem.

Since the same range is used for errors for both continuous and discrete adjoints, the magnitudes of the errors may be compared easily and it is obvious that for all cases, the magnitude of the error in the corrected functional based on the continuous adjoint is less than the corresponding value for the discrete adjoint. This is expected behavior as a result of solving the continuous adjoint PDE. Furthermore, the convergence rates are faster for all cases using the continuous adjoint. In conclusion, for advection, significant improvement is obtained by correcting the functional by the truncation error based on the smoothed solution, whether the continuous or discrete adjoint is used, although the continuous adjoint is more successful.

6.1.2. Poisson equation Consider the 2D Poisson equation:

u = f in D u = bD on ∂ DD nˆ · ∇ u = bN on ∂ DN

(26)

where bD is the Dirichlet boundary condition on ∂ DD and bN is the Neumann boundary condition on ∂ DN , with ∂ DD ∪ ∂ DN = ∂ D and ∂ DD ∩ ∂ DN = 0. Similar to the advection problem, the left hand side of Eq. (26) is multiplied by z, integrated over the domain D and integrated by parts to obtain the continuous adjoint equation:

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

413

only quadratic. This is the behavior of diffusive fluxes for thirdorder finite volume unstructured schemes and as the third-order functional is not third-order, as shown in Fig. 8, correcting it with fourth-order flux integration is not helpful either. Summarizing, for Poisson, significant improvement is obtained by correcting the functional by the truncation error based on the smoothed solution, when either the continuous or discrete adjoint is used, particularly for C s J24 , although the improvement is more significant when the continuous adjoint is used. 6.2. Systems of equations Our results for model problems show that the smoothing scheme is successful in correcting the functionals for both convective and diffusive fluxes which suggests that it can be extended to real problems involving system of equations. Both inviscid and viscous flows are considered in this section. Since correcting the functional by the truncation error based on the non-smoothed solution is not helpful for advection and Poisson, we will not show the results for the non-smoothed solution anymore. 6.2.1. Euler equations Consider the linearized steady-state Euler equations: Fig. 7. Exact solution of the Poisson problem.

    (u, z )D + u, nˆ · ∇ z ∂ D + nˆ · ∇ u, −z ∂ D D    N = (u, z )D + nˆ · ∇ u, z + u, −nˆ · ∇ z ∂D ∂D D

LU =

Cu = nˆ · ∇ u

Bu = nˆ · ∇ u,

⎛

⎞ ρ uy ⎜ρ u2x + p⎟ c ⎜ ρ ux uy ⎟ ⎟ ⎜ ⎟ Fxc = ⎜ ⎝ ρ ux uy ⎠, Fy = ⎝ρ u2y + p⎠ ρ ux H ρ uy H

N

on ∂ DN

C ∗ z = nˆ · ∇ z

B z = −nˆ · ∇ z, ∗

∗

on ∂ DD

C z = −z



on ∂ DN

for the continuous adjoint problem. The output functional is defined as:

ˆ

J=

D

ˆ

uzdA +

ˆ =

D

ˆ zudA +



∂ DD ∂ DD



u nˆ · ∇ z ds +





nˆ · ∇ u zds +

ˆ

ˆ



∂ DN ∂ DN





nˆ · ∇ u (−z )ds



−nˆ · ∇ z uds

(27)

The primal problem is defined on a (1 × 1) square domain as

∂ 2u ∂ 2u π3 + = − sin (π x ) sin (π y ), u = sin (π x ) on y = 1 4 ∂ x2 ∂ y2 ⇒ uexact =

π 8

sin (π x ) sin (π y ) +

ρ ux

⎞

⎛

(30)

where U = (ρ ρ ux ρ uy ρ E )T is the conserved solution vector. Eq. (29) is multiplied by Z, integrated over the domain D and integrated by parts to obtain the continuous adjoint equation:

for the primal problem and L∗ z = z in D and

B∗ z = z,

∂Fc

∂Fc

on ∂ DD

Cu = u

(29)

where Ax = ∂ Ux and Ay = ∂ Uy in which Fx and Fy are the convective flux functions:

Comparing to Eq. (12), we end up with Lu = u in D and

Bu = u,

∂ ∂ ( A U ) + ( Ay U ) = f ∂x x ∂y

1 sin (π x ) sinh (π y ) sinh (π ) (28)

and the exact solution is shown in Fig. 7. Setting the adjoint source 5 term as g = π2 x(1 − x )y(1 − y ), the adjoint boundary condition is obtained as h = x(1 − x ) on y = 0 and zero elsewhere by Eq. (12). The output functional is defined as:

J = (u, g)D + (Cu, h )∂ DD = (u, g)D + (nˆ  ∇ u, h )∂ DD We correct the functional using the adjoint solutions; the results are shown in Fig. 8. Similar to the advection problem, correcting the functional by the non-smoothed truncation error is not helpful and on the other hand, the corrected functional by smoothed truncation error and the continuous adjoint solution provides the best convergence rate. For the Poisson problem, we are comparing the error to the fourth-order functional instead of the thirdorder one since the convergence of the third-order functional is

 ∂ ∂ A U + A U , Z ( ) ( ) ∂x x ∂y y  D ∂ Z ∂ Z = U, −ATx − ATy + ( (nx Ax + ny Ay )U, Z )∂ D ∂x ∂y D

(31)

Comparing to Eq. (12), we end up with L∗ Z = −ATx ∂∂ Zx − ATy ∂∂ Zy for the continuous adjoint problem. The boundary operators are different for supersonic and subsonic inflow and outflow. Details of different boundary conditions and operators are described by Pierce and Giles [12] and Alauzet and Pironneau [3]. Since lift and drag are the most common functionals of interest in several cases, the boundary operators corresponding to these functionals on walls are illustrated here. The output functional corresponding to drag or lift is set as:

ˆ

J=

∂ Dwall

 ds pnˆ · ψ

(32)

where ψ = (cos α , sin α )T for drag and ψ = (− sin α , cos α )T for lift, and α is the angle of attack. Using the primal boundary operator in the right hand side of Eq. (31) for the slip wall boundary  · nˆ = ux nx + uy ny = 0 and multiplying by the adjoint socondition u lution,

( (nx Ax + ny Ay )U, Z )∂ D =

 

0

z1

pnx z2

pny z3

z4

T

0 ,

T

= pnx z2 + pny z3

Comparing to the functional of interest, Eq. (32), we obtain the continuous adjoint boundary condition as

 nx z2 + ny z3 = nˆ · ψ

(33)

414

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

Fig. 8. Convergence history for Poisson.

Supersonic flow The supersonic vortex is an Euler problem with an exact solution for which the exact functional value can therefore be calculated. The exact solution is

    γ 1−1 R2 γ −1 ρ = 1+ 1 − i2 Mi2 2

R

yRi Mi ux = R2 −xRi Mi uy = R2

ργ γ

p=

where we take Ri = 2.0 as the inner radius and Mi = 2.0 is the inlet Mach number. Considering drag or lift on the inner wall of the supersonic vortex and using Eq. (32) for calculating the functional based on the exact pressure distribution, the exact value for the functional is obtained as:

ˆ

J= =

inner wall

1

γ

ˆ 0

π 2



 ds pnˆ · ψ

 1 + 2 (γ − 1 ) 1 −

R2i R2

 γ γ−1

 dθ = 1.42857 Ri nˆ · ψ (34)

It should be noted that the same magnitude is obtained for both drag and lift. The adjoint solution distributions based on drag and lift functionals on the lower wall are depicted in Fig. 9 for both continuous and discrete adjoint solutions. Similar to the advection test case, the continuous adjoint solution distribution is smoother than the discrete adjoint. While both functionals show high sensitivities near the wall, the drag functional adjoint shows high sensitivities along a longer extent than the lift functional due to the need to accurately convect information to the outflow region, where a large amount of the drag force is produced. For the same solution component, x-momentum, the convergence results are shown in Fig. 10 comparing the discretization error for second order accurate solution and the exact truncation

error. As mentioned earlier, the magnitude of the truncation error for the unstructured mesh is larger than the discretization error and the convergence rate is slower. The estimation of the truncation error using the p-TE method is also compared to the smoothed truncation error and the truncation error magnitude is smaller for the smoothed estimate. Functional correction based on these adjoint solutions is applied and the results are shown in Fig. 11. The corrected functional based on the smoothed solution, C s J23 , has smaller error compared to the second-order functional and it converges cubically. Fourthorder convergence is obtained for C s J24 based on the continuous adjoint for both functionals; however, using the discrete adjoint for C s J24 , lift and drag have different behaviors. Although the convergence for drag is even better than quartic, the convergence rate for lift is only 2.54. Generally, using the continuous adjoint for functional corrections provides us with more consistent results, as it is less sensitive to the errors arising from discretization. Correction based on third order solution has the same behavior and the best convergence rate is based on the continuous adjoint solution. Consequently, as a general conclusion, correcting the functional based on the non-smoothed solution is not helpful at all and smoothing the solution and computing the correction term based on the smoothed solution is capable of improving the functional significantly. More accurate results with the expected convergence behavior are obtained when the continuous adjoint solution is used for correction. Subsonic flow The capability of our smoothing scheme in correcting the functional for an Euler problem with an exact solution was shown for the supersonic vortex. In this section, output functionals of drag and lift are calculated and corrected for the subsonic inviscid flow around NACA-0012 with Ma∞ = 0.5 and angle of attack α = 2◦ . Since this problem does not have an exact solution, the solution on a very fine mesh is considered as the exact solution. The capability of our solver in approximating the drag and lift for inviscid flow around NACA-0012 with Ma∞ = 0.5 and α = 2◦ was shown in the first International Higher-Order Workshop [1] and the secondorder solution on a mesh with characteristic cell size nine times

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

415

Fig. 9. x−momentum adjoint solutions for supersonic vortex.

smaller in area than the finest mesh for which correction is applied is considered as the exact solution to the problem. The density component of the solution is depicted in Fig. 12 and the farfield boundary is a circle of radius 500 chords centered at the leading edge of the airfoil and consequently, the errors arising from boundary placement can be neglected. Functional correction based on these adjoint solutions is applied and the results are summarized in Fig. 13. Similar to the other cases being tested so far, the corrected functional using the continuous adjoint solution provides the best convergence rate and the magnitude of the error is much smaller than the third order solution for both lift and drag functionals.

6.2.2. Navier–Stokes equations Consider the linearized steady-state Navier–Stokes equations:

  ∂ ∂U ∂U AxU − Dxx − Dxy ∂x ∂x ∂y   ∂ ∂U ∂U + AyU − Dyx − Dyy = f ∂y ∂x ∂y

LU =

(35)

where Ax and Ay are the derivatives of the convective flux function, Eq. (30), with respect to the solution as defined before and the other Dij matrices are the derivatives of viscous flux with respect

416

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

Fig. 10. Convergence results for reconstruction scheme, supersonic vortex.

to the solution gradients

∂Fv ∂ Fxv  ∂U  , Dyx =  ∂yU  ∂ ∂x ∂ ∂x v ∂ Fyv ∂F Dxy =  ∂xU  , Dyy =  ∂ U  ∂ ∂y ∂ ∂y

6.3. Remarks on computational cost

Dxx =

and Fv represents the viscous flux function. Eq. (35) is multiplied by Z, integrated over the domain D and integrated by parts to obtain the continuous adjoint equation:

∂Z ∂Z − ATy ∂x ∂y     ∂ ∂Z ∂Z ∂ ∂Z ∂Z − DTxx + DTyx − DTxy + DTyy ∂x ∂x ∂y ∂y ∂x ∂y

L∗ Z = −ATx

for the continuous adjoint problem. Inflow and outflow boundary conditions and operators are defined by Pierce and Giles [12]. Since there is not any commonly-used Navier–Stokes problem with an exact solution in the literature, we use a manufactured solution [28] to verify the performance of our smoothing scheme for correcting the functional. We use a manufactured solution for both the primal [32] and adjoint problems. The general form of the manufactured primal solution is

U (x, y ) = a0 + ax sin (bx π x + cx π ) + ay sin (by π y + cy π ) +axy sin (bxy π y + cxy π )

(36)

and the general form of the manufactured adjoint solution is

Z (x, y ) = a0 [sin (π (1 − x ) ) sin (π (1 − y ) ) sin (bx π x + cx π ) sin

(by π y + cy π ) sin (bxy π y + cxy π ) + (1 − x )(1 − y )] (37) where the coefficients for four components are as illustrated in Table 1. The inlet Mach number is set as Mai = 2.5 and the Reynolds number is defined as

Re =

the manufactured primal and adjoint solutions are calculated easily and the corresponding weight boundary operators are found by integration by parts, similar to the Dirichlet and Neumann boundary operators obtained for the Poisson problem. Functional correction based on these adjoint solutions is applied and the results are shown in Fig. 15. Similar to the conclusion for the previous test case, using the continuous adjoint for correction improves the functional more and provides us a faster convergence rate as well. For this test case, the manufactured solution is used for both primal and adjoint problems, and consequently calculating the error of the adjoint solution is easy. Fig. 16a depicts the error for both discrete and continuous solutions. Not surprisingly, the magnitude of the error in calculating the continuous adjoint is less than the discrete adjoint and the difference becomes smaller by decreasing the mesh size and both converge quadratically. In the beginning of this section, verifying adjoint consistency was deferred for consideration with the Navier–Stokes equations, because this is a system of equations with both convective and diffusive terms. Adjoint consistency for the discrete adjoint solution for this test case is depicted in Fig. 16 which shows that almost the same functional value is obtained either based on the primal solution or the adjoint solution and the difference decreases with mesh refinement, implying that the discretization is adjoint consistent.

ρi u i L = 800 μi

where the characteristic length is L = 1.0, ρi = 1.0 and μi = 1.0. The y−velocity of the primal solution is shown in Fig. 14. Since we are using a manufactured solution for both primal and adjoint problems, the functional is directly obtained by the adjoint definition, Eq. (12). The source term and boundary conditions of

Output functional correction for unstructured meshes converges to the exact value with the expected convergence rate, similar to the structured mesh cases, but by estimating the truncation error based on a C1 interpolated smoothed solution, higher order convergence of the functional is obtained by using the lower order solution. Hence, paying the extra cost of computing the higher-order solution to obtain a more accurate functional which converges to the exact value with higher-order rate is not required and the higher-order convergence is obtained. Smoothing the solution using C1 interpolation, however, is an extra computational cost compared to the lower-order solution. The extra computational cost includes solving the (21 × 21) system of equations corresponding to 21 constants required for 21 equations of Argyris element. It should be noted that this system of equations is solved only once and all elements are mapped into a single Argyris reference element and accordingly this post-processing step is cheap in terms of time and memory. A comparison between the CPU time required for second, third and fourth-order solutions and finding the C1 interpolation of the solution is shown in Table 2. Two different problems are considered: a scalar Poisson problem and an Euler system of equations. These two problems are described earlier in this section. Investigation of Table 2 shows that the post processing time required for C1 interpolation of the solution is small compared to the solution time, and scales linearly with problem size, whereas the solution time grows faster than linearly. The CPU time required for C1 interpolation of the third order solution is more than interpolation of the second order one as finding the solution and first and second derivatives of the solution is based on reconstruction and third order reconstruction needs more time compared to second order. There is the same behavior for fourth order compared to third order solution. Another point that needs to be considered about computational cost is the time comparison between continuous and discrete adjoint solutions. For the discrete adjoint, the computational cost includes finding the transpose of the Jacobian and the derivative of the functional with respect to the solution and then solving the system of equations once for which the computational cost is of the order of one linear system solution. On the other hand, for the continuous adjoint, the computational time is at the same or-

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

417

Fig. 11. Convergence history for supersonic vortex. Table 1 Coefficients for manufactured primal and adjoint solution used in Eq. (36) and (37). c is the speed of sound in SI units and is used for non-dimensionalization .

ρ ux uy p

a0

ax

bx

cx

ay

by

cy

axy

bxy

cxy

1.0 800/c 800/c 105 /c2

0.15 50/c −75/c −2 × 104 /c2

2.0 0.66 1.66 0.5

0.33 0.5 −0.5 0.5

−0.1 −30/c 40/c 5 × 104 /c2

1.0 1.5 3.0 1.0

−0.2 −0.125 −0.0 −0.5

−0.05 −60/c 60/c −1 × 104 /c2

1.0 1.5 1.5 0.75

0.25 0.125 0.125 −0.25

der of computational time required for the primal problem. Hence, the total time for correcting the second-order functional using our method with the continuous solution adjoint is about twice as large as solving the primal problem, while the computational time for correction using the discrete adjoint is the equal to the time for solving the primal problem plus one non-linear iteration.

7. Concluding remarks It is known that the truncation error on unstructured meshes has larger magnitude compared to structured meshes and the rough modes in the truncation error are dominant. In this paper, we have presented a C1 interpolation scheme to obtain continu-

418

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

Table 2 CPU time (sec) comparison for solution and C1 interpolation of the solution. 2nd order

3rd order

#CV s

solution

solution

C1 interpolation

solution

C1 interpolation

1260 5004 19862 80430

Poisson 0.143 0.795 9.305 76.875

0.353 1.635 14.418 130.210

0.069 0.306 1.257 5.092

0.615 2.635 18.495 158.069

0.201 0.852 3.251 13.192

1044 4408 16256 65502

Euler 0.197 0.959 7.787 94.411

0.448 2.072 10.993 202.422

0.063 0.291 1.111 4.446

0.843 3.915 18.658 264.728

0.171 0.723 2.718 11.626

Fig. 12. x−momentum distribution for NACA-0012, Ma∞ = 0.5, α = 2◦ .

Fig. 13. Convergence history for NACA-0012, Ma = 0.5, α = 2◦ .

4th order

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

Fig. 14. Exact solution of Navier–Stokes problem, y−velocity, Re = 800 , Ma = 0.5 .

ous solution and gradients at the quadrature points on the faces of neighboring cells. By using the continuous interpolation, the problem of having discontinuous jumps of the coefficients of the leading terms in the Taylor series expansion of the error from one control volume to another is solved and we can obtain a single value for the flux vector from two adjacent cells. Since we are using the p−truncation error approach for estimating the truncation error, which computes the higher-order flux integration based on the lower-order solution, the smooth solution and consequently the smooth flux vector is able to provide us a smoother estimate of the truncation error. The truncation error estimate when multiplied by the adjoint solution is used to compute a correction term for an output functional of interest. We combine the smooth truncation error estimate with both continuous and discrete adjoints and compare their performance.

419

Our results show that using a truncation error estimate based on the non-smooth solution does not improve the value of the functional, neither in terms of the convergence rate nor in terms of the error magnitude, for all problems tested. Our test cases are model problems with convective and diffusive fluxes, the advection and Poisson equations, respectively, and real flows for both inviscid and viscous flows. Although using the truncation error based on the non-smoothed solution is not helpful, using the C1 interpolated solution improves the functional significantly. The second-order solution is used to find a continuous flux integral based on the third or fourth-order reconstruction giving third or fourth-order convergence for the functional while only solving the primal and adjoint problems at second-order. Using fourth-order flux integration based on the third-order solution provides us with similar results and for both cases, the magnitude of the functional error is also less than the magnitude of the uncorrected functional. More accurate corrected functionals are obtained when correction is based on the continuous adjoint since a more accurate PDE is solved compared to solving a linear system which gives the discrete adjoint. As a result more consistent and accurate results for different problems with different orders are obtained using the continuous adjoint. However, a more accurate solution is compensated by more computational costs for the continuous adjoint. Computational time for the continuous adjoint is about the same as computational time for the primal problem. For the discrete adjoint, the computational time is the time of solving a linear system, so about the cost of one time step for the primal. This significant improvement in functional correction is obtained just by smoothing the primal solution and finding the truncation error based on this smoothed solution and then multiplying the truncation error by the adjoint solution of the same order as the primal solution. It is not necessary to smooth the adjoint solution to obtain significant improvement in convergence rates. The smoothing scheme we developed can be extended to any other two-dimensional problem without any difficulties as it worked properly for both viscous and inviscid flows which involve diffusive and convective terms. Extending that to 3D problems involves finding the proper polynomials for the 3D reference element, similar to the Argyris element for 2D, and then solving the

Fig. 15. Convergence history for Navier–Stokes, Re = 800 , Ma = 0.5.

420

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

Fig. 16. Adjoint solution for Navier–Stokes. Table 3 Non-zero coefficients in Eq. (21) for the Argyris reference element. c

ϕ0 ϕ1 ϕ2 ϕ3 ϕ4 ϕ5 ϕ6 ϕ7 ϕ8 ϕ9 ϕ 10 ϕ 11 ϕ 12 ϕ 13 ϕ 14 ϕ 15 ϕ 16 ϕ 17 ϕ 18 ϕ 19 ϕ 20

x

y

x2

xy

y2

x3

x2 y

−10 −6

1 1

xy2

x4

−10 −11 −11

1

y3

18 1.5

−4

1

x2 y2

15 8

−6

−1.5

0.5

x3 y

−4

5 −1.5

0.5

−15 7

10 −4 −5

14 −1

0.5

−3

1 10 −5 −4 1 0.5

−32

16 −16

system of equations arising from the basis functions once to obtain the constants. Acknowledgment This work was supported by ANSYS Canada and the Canadian Natural Sciences and Engineering Research Council under Collaborative Research and Development Grant 431183-12. Appendix A. C1 Interpolation of solution The reference Argyris element is a right isosceles triangle and the vertices are set as shown in Fig. 3. As mentioned, 21 basis functions, ϕ i , are required. The basis functions are defined as Eq. (22) which implies that 21 constants are required for each. The

−30 10 10 −1.5 10 −1.5 15 −3.5 18.5 0.25 −3.5 1.25 15 18.5 −3.5 1.25 −3.5 0.25 −32 √ 8 2 32

xy3

y4

x5

15 18

x4 y

−6 −3 −8

8 −0.5

−2

5 1.5 6 −3

−8 0.5 2 −15 14 7 −3 −1 16 32

x3 y2

x2 y3

30 1 −10 1.5 −6 1 −15 3.5 −18.5 −0.25 3.5 −0.75 −15 −13.5 3.5 −1.25 2.5 −0.25 32 √ −8 2 −16

30 −10 1 1 −6 1.5 −15 3.5 −13.5 −0.25 2.5 −1.25 −15 −18.5 3.5 −0.75 3.5 −0.25 16 √ −8 2 −32

xy4

−6 −8 −3 −2 −0.5

6 −8 −3 2 0.5

−16

non-zero constraint on the basis functions are:

ϕ0 |(0,0) = 1, ∂ϕ1  = 1, ∂ x  ( 0,0 )  ∂ϕ2  = 1, ∂ y  ( 0,0 )  ∂ 2 ϕ3  = 1, ∂ x2  ( 0,0 )  ∂ 2 ϕ4  = 1, ∂ x∂ y  ( 0 , 0 )  ∂ 2 ϕ5  = 1, ∂ y2  ( 0,0 )

ϕ6 |(1,0) = 1, ∂ϕ7  = 1, ∂ x  ( 1,0 )  ∂ϕ8  = 1, ∂ y  ( 1,0 )  ∂ 2 ϕ9  = 1, ∂ x2  ( 1,0 )  ∂ 2 ϕ10  = 1, ∂ x∂ y  ( 1 , 0 )  ∂ 2 ϕ11  = 1, ∂ y2  ( 1,0 )

y5

ϕ12 |(0,1) = 1 ∂ϕ13  =1 ∂ x  ( 0,1 )  ∂ϕ14  =1 ∂ y  ( 0,1 )  ∂ 2 ϕ15  =1 ∂ x2  ( 0,1 )  ∂ 2 ϕ16  =1 ∂ x∂ y  ( 0 , 1 )  ∂ 2 ϕ17  =1 ∂ y2  ( 0,1 )

M. Sharbatdar et al. / Computers and Fluids 140 (2016) 406–421

   ∂ϕ18 ∂ϕ18  ∂ϕ18 nξ + nη  = − =1 ∂x ∂y ∂ y  1 ,0 ( 12 ,0 ) (2 )    √   ∂ϕ19 2 ∂ϕ19 ∂ϕ19  ∂ϕ19 nξ + nη  = + =1 ∂x ∂y 2 ∂x ∂y  1 , 1 ( 12 , 12 ) (2 2)      ∂ϕ20 ∂ϕ ∂ϕ20 20  nξ + nη  = − =1 ∂x ∂y ∂ x  0, 1 (0, 12 ) ( 2) 

This system of equations is solved once and the corresponding coefficients for the basis functions are as listed in Table 3. References [1] First International Workshop on High-Order CFD Methods. 2012. http://people. ku.edu/∼z651w035/hiocfd.html; Sponsored by the AIAA Fluid Dynamics Technical Committee, the US Air Force Office of Scientific Research, and the Deutsches Zentrum fur Luft- und Raumfahrt. [2] Ainsworth M, Oden JT. A posteriori error estimation in finite element analysis. Wiley-Interscience, John Wiley and Sons; 20 0 0. [3] Alauzet F, Pironneau O. Continuous and discrete adjoints to the Euler equations for fluids. Int J Numer Methods Fluids 2012;70:135–57. [4] Argyris JH, Fried I, Scharpf DW. The tuba family of plate elements for the matrix displacement method.. Aeronaut J R Aeronaut Soc 1968;72:701–9. [5] Becker R, Rannacher R. Weighted a posteriori error control in finite element methods. In: Proceedings of ENUMATH-97, Heidelberg; 1998. [6] Ciarlet PG. The finite element method for elliptic problems. Soc Ind Appl Math (SIAM) 2002. [7] Diskin B, Thomas JL. Accuracy analysis for mixed-element finite-volume discretization schemes. Tech. Rep.. National Institute of Aerospace; 2007. [8] Diskin B, Thomas J. Notes on accuracy of finite-volume discretization schemes on irregular grids. Appl Numer Math 2010;60:224–6. Rebuttal of Svard 2008 [9] Diskin B, Thomas JL, Nielsen EJ, Nishikawa H, White JA. Comparison of node– centered and cell-centered unstructured finite-volume discretizations: viscous fluxes. AIAA J 2010;48(7):1326–38. [10] Fidkowski K, Darmofal D. Reviw of output-based error estimation and mesh adaptation in computational fluid dynamics. AIAA J 2011;49(4):673–94. [11] Giles MB, Pierce NA. An introduction to the adjoint approach to design. Flow, Turbul, Combust 20 0 0;65(3–4):393–415. [12] Giles MB, Pierce NA. Adjoint Equations in CFD: Duality, Boundary Conditions and Solution Behaviour. In: 13th Computational Fluid Dynamics Conference; 1997. [13] Jalali A, Sharbatdar M, AOllivier-Gooch C. Accuracy analysis of unstructured finite volume discretization schemes for diffusive fluxes. Comput Fluids 2014;101:220–32. [14] Jalali A, Ollivier-Gooch C. Accuracy assessment of finite volume discretizations of diffusive fluxes on unstructured meshes. In: Proceedings of the Fiftieth AIAA Aerospace Sciences Meeting; 2012. AIAA Paper 2012-0608 [15] Jameson A, Martinelli L, Pierce NA. Optimum aerodynamic design using the Navier–Stokes equations. Theor Comput Fluid Dyn 1998;10(1–4):213–37. [16] Jameson A. Aerodynamic design via control theory. J Sci Comput 1988;3(3):233–60. [17] Jameson A, Shankaran S, Martinelli L. Continuous adjoint methods for unstructured grids. AIAA J 2008;46:1226–39. [18] Kim H, Nakahashi K. Unstructured Adjoint Method for Navier–Stokes Equations. JSME Int J Series B 2005;48:202–7.

421

[19] Kim H, Sasaki D, Obayashi S, Nakahashi K. Aerodynamic optimization of supersonic transport wing using unstructured adjoint method. AIAA J 2001;39:1011–20. [20] Larson MG, Barth TJ. A posteriori Error Estimation for Discontinuous Galerkin Approximations of Hyperbolic Systems. NASA Technical Report; 1999. [21] Leicht T, Hartmann R. Anisotropic mesh refinement for discontinuous galerkin methods in two-dimensional aerodynamic flow simulations. Int J Numer Methods Fluids 2008;56:2111–38. [22] Michalak C, Ollivier-Gooch C. Globalized matrix-explicit newton-GMRES for the high-order accurate solution of the euler equations. Comput Fluids 2010;39:1156–67. [23] Mohammadi B, Pironneau O. Shape optimization in fluid mechanics. Annu Rev Fluid Mech 2004;36:11.1–11.25. [24] Nejat A. A higher-order accurate unstructured finite volume Newton-Krylov algorithm for inviscid compressible flows. The University of British Columbia, Department of Mechanical Engineering; 2007. [25] Nemec M, Aftosmis MJ. Adjoint error estimation and adaptive refinement for embedded-boundary cartesian meshes. In: Proceedings of the Eighteenth AIAA Computational Fluid Dynamics Conference; 2007. [26] Nemec M, Aftosmis MJ, Wintzer M. Adjoint-based adaptive mesh refinement for complex geometries. Forty-sixth AIAA Aerospace Sciences Meeting; 2008. AIAA 2008-725. [27] Nishikawa H. Beyond interface gradient: A general principle for constructing diffusion scheme. In: Proceedings of the Fortieth AIAA Fluid Dynamics Conference and Exhibit; 2010. AIAA paper 2010–5093. [28] Oberkampf WL, Roy CJ. Verification and validation in scientific computing.. Cambridge University Press 2010. [29] Ollivier-Gooch C, Nejat A, Michalak C. On obtaining and verifying high-order finite-volume solutions to the euler equations on unstructured meshes. AIAA J 2009;47(9):2105–20. [30] Ollivier-Gooch C, Roy C. Reducing truncation error on unstructured meshes by vertex movement. In: Proceedings of the Forty-Second AIAA Fluid Dynamics Conference; 2012. [31] Ollivier-Gooch CF, Van Altena M. A high-order accurate unstructured mesh finite-volume scheme for the advection-diffusion equation. J Comput Phys 2002;181(2):729–52. [32] Phillips T, Roy C, Borggaard J. Error Transport Equation Boundary Conditions for the Euler and Navier-Stokes Equations. 52nd Aerospace Sciences Meeting, AIAA SciTech; 2014. [33] Pierce N, Giles M. Adjoint recovery of superconvergent functionals from PDE approximations. SIAM Rev 20 0 0;42(2):247–64. [34] Pierce N, Giles M. Adjoint and defect error bounding and correction for functional estimates. J Comput Phys 20 04;20 0(2):769–94. [35] Pironneau O. On optimum design in fluid mechanics.. J Fluid Mech 1974;64:97–110. [36] Qiao Z, Yang X. Wing Design by Solving Adjoint Equations. 40th AIAA Aerospace Sciences Meeting & Exhibit; 2002. [37] Roache PJ. Code verification by the method of manufactured solutions. J Fluids Eng 2002;124(1):4–10. [38] Roe PL. Characteristic-based schemes for the euler equations. In: Annual Review of Fluid Mechanics, 18. Annuals Reviews, Inc.; 1986. p. 337–65. [39] Sharbatdar M, Ollivier-Gooch C. Eigenanalysis of truncation and discretization error on unstructured meshes. In: 21st AIAA Computational Fluid Dynamics Conference; 2013. [40] Venditti DA, Darmofal DL. Adjoint error estimation and grid adaptation for functional outputs: application to quasi-one-dimensional flow. J Comput Phys 20 0 0;164(1):204–27. [41] Venditti DA, Darmofal DL. Grid adaptation for functional outputs: application to two-dimensional inviscid flows. J Comput Phys 2002;175(1):40–69.