Introduction

The approximation error is omnipresent at the numerical solutions of Partial Differential Equations (PDE) due to the discretization at numerical statements. The error estimation is of the high current interest in view of the need for the verification of software and numerical solutions. For example, the corresponding issues are stated in Computational Fluid Dynamics (CFD) in the form of standards [1,2]. Let's discuss main approaches to the estimation of the approximation error. We consider a PDE system in the operator form

f u A = ~ (1

) and corresponding discrete operatorh h h f u A =(2)

that approximates the system on some grid.

In further presentation we consider the numerical solution uh to be a grid function (vector u h ∈RM, M is the number of grid points), M h R u ∈ ~ to be the projection of true solution onto the grid,

h h h u u u − = ∆

to be the true approximation error, ∆uh to be some estimate of this error. L2 -based norm is used for a comparison of these vectors. We may also use the set of numerical solutions

M i h R u ∈ ) (

obtained by independent numerical algorithms (i=1…K) is the number of used algorithm). Two main options to the estimation of the approximation error exist.

A priori error estimation n Ch u ≤ ∆ (3) is commonly used at the design and the theoretical analysis for the determination of the convergence order. Herein h is the step of discretization, n is the order of approximation, C is an unknown constant. This estimation is usually related to the truncation error (source term of differential approximation [3]) that may be formally expressed as

1 1+ + ∞ = ∂ ∂ = ∑ m m m m n m x u h C u δ (4)

for PDE systems of the first order. This expression contains the infinite number of terms so, the first (minor order) term is commonly used. It may be computed by many ways including the special postprocessor [4].

A priori error estimates have an universal form. Unfortunately, the unknown constant prevents it from to be used in applications.

A posteriori error estimator usually has the form

) ( h u E u ≤ ∆ (5)

and is determined by some computable error indicator E(uh). This estimator depends on the previously computed numerical solution uh and, thus, has a minor generality. Fortunately, it may be applied to practical computations since has no unknown constants.

The highly efficient technique is developed for a posteriori error estimation in the domain of the finiteelement analysis [5,7,8]. In accordance with [5], the quality of a posteriori error estimation may be expressed via the effectivity index of estimator that is equal to the relation of the estimated error norm to the true error norm:

2 2 ) ( ) ( ) ( ~L i h L i h j eff u u I ∆ ∆ = (6)

One may treat the norms of the true error and estimation error as the radii of hyperspheres

I

means that the hypersphere, containing the true error, belongs to the hypersphere defined by the estimator. So, in order to provide the reliable estimation, the effectivity index should be greater the unit. On the other hand, the estimation should be not too pessimistic, so the value of the effectivity index should be not too great. For the finite element methods, used in the domain of elliptic equations (usually engendering highly regular solutions), the acceptable range, according [5], is 3 1

) ( ≤ ≤ i eff I .

However, the boundaries of this inequality seem to be dependent on the problem at hand. Numerical tests for CFD domain demonstrate the efficiency index to belong the range (0.3,5). For the nonlinear problems containing discontinuities (that is common case for CFD) the progress of a posteriori error estimation in the rigorous form of inequality ( 5) is limited.

As an alternative, some less rigourous methods are employed. These methods provide the estimation of ∆uh without any strict inequality. Significant number of estimators do not met the condition 1

)

( ≥ i eff I .

We note such error estimators as nonstrict (weak) ones.

The first domain of these approaches forms the defect correction methods [9][10][11]. Some part of these methods [9,10] apply some approximation of the truncation error δuh in order to disturb the main system. The additional equation for the error transformation occurs

h h h h u u A δ = ∆ ) ( .(7)

These methods are rather laborious since imply coding, debugging and solving of an additional problem.

A bit less laborious version of defect correction methods [11] have the appearance

h h refined h h h f u u A δ + = ) ( ,(8)

that imply the disturbing of the main problem by the source term, which approximates the truncation error.

Another branch of nonstrict a posteriori error estimation methods has a non-intrusive form of certain postprocessor that significantly reduces efforts for coding and debugging. It may be based on the Runge rule [5], Richardson extrapolation (RE) [13,14], Inverse Problem based approach (IP) [15] or ensemble based methods (EM) [16][17][18][19].

The heuristic rule by C. Runge [5] is the basis of commonly used stopping criterion by merging of solutions at the mesh refinement.

The standards for verification and validation [1,2] recommend the Richardson extrapolation (RE) as the main tool for the verification. Richardson extrapolation provides the pointwise approximation of the error field, unfortunately, at the cost of the high computational burden [13,14]. RE provides some generalization of the Runge's rule.

The Inverse Problem based approach (IP) [15] enables the pointwise information on the error.

The computationally cheap approach to a posteriori error estimation that is based on the ensemble of numerical solutions obtained by independent methods is offered by [16][17][18][19]. However, these approaches do not provide the pointwise information on the approximation error.

Thus, the computationally inexpensive nonstrict a posteriori estimation of approximation error is of the major interest in CFD from the viewpoint of verification of codes and solutions. The simultaneous use of several nonstrict methods may have some prospects from the viewpoint of reliability increasing.

Runge rule

From the historic viewpoint the first a posteriori error indicator is based on the heuristic rule by C. Runge [5]. If the difference between two approximate solutions computed on a coarse mesh uh and the refined mesh uh,ref is small, then both are assumed to be close to the exact solution. The Runge's rule can be considered as the first

a posteriori error indicator |ε(uh)-ε(uh,ref)|=ERunge(uh-+ + + = − = ∆ j j j h h h C h C h C u u u (9)

where j1, j2, j3, … are positive (sometimes noninteger) numbers ordered in accordance with the magnitude (for example, [12]).

The accuracy for the error estimation by Runge's rule has the lowest order ) ( 1 j h O and remains unresolved.

Richardson extrapolation

Richardson extrapolation is based on the first term of expansion (9) n q q Ch u u + = ) ( and operates if the asymptotic range is achieved for several grids hq (Ck, n are assumed to be constant that should be verified numerically by expanding the set of grids) [13,14]. For CFD problems containing shock waves and contact lines [12] the error order is not constant over the flowfield and depends on the type of flow structure. So, one should to extend RE for estimation of the local order of convergence.

The pointwise (m is the grid point number) results of numerical computation for three meshes of different steps hq may be presented as:

m n k m m h C u u 1 ) 1 ( ~+ = m n k m m h C u u 2 ) 2 ( ~+ = (10) m n k m m h C u u 3 ) 3 ( ~+ = .

This system is defined for the most rough grid and may be resolved regarding

m m m n C u , , ~

by several methods described in [13,14]. The relations (10) are valid, if Cm is independent on h and higher order terms may be neglected, that is, the solution is in the asymptotic range. This statement implies at least four consequently refined grids. If the asymptotic range is not confirmed on these grids, the additional refinement is necessary.

So, the Richardson extrapolation provides the pointwise approximation for the error field at the cost of the high computational burden.

The accuracy of RE for the error estimation has the appearance ) ( 2 j h O and remains unresolved quantitatively that excludes estimates by inequalities of the type (5).

Approximation error estimation using Inverse Problems

The nonstrict (weak) form of the approximation error estimation by the Richardson extrapolation uses the set of numerical solutions obtained by the same algorithm for consequently refined grids. On other hand, one may consider an ensemble of numerical solutions

) ( , ) ( ~i m m h i m u u u ∆ + =

, obtained by K independent algorithms (of different structure, for example, of different approximation order) on the same grid. The projection of an exact (unknown) solution on the grid point m is denoted as m h u , ~, the approximation error (also unknown) for i-th solution is denoted as

) (i m u ∆ .

The differences of solutions are equal to the differences of the approximation errors and, hence, contain some information regarding the unknown errors

) (i m u ∆ : ) ( ) ( ) ( ) ( , j m i m j m i m m ij u u u u d ∆ − ∆ = − =(11)

We treat this information in accordance with approach described by [15] in order to determine the approximation error

) (i m u ∆ . One may obtain 2 / ) 1 ( − ⋅ = n n N

relations on the set of n numerical solutions:

m i j m ij f u A , ) ( = ∆ . (12)

The summation over the repeating indexes is used elsewhere from this point. For the minimum set of data (three solution) equation ( 12) assumes the form

          − − − =           =           ∆ ∆ ∆           − − − ) 3 ( ) 2 ( ) 3 ( ) 1 ( ) 2 ( )1 ( , 3 , 2 , 1 ) 3 (+ ∆ = ∆ ) ( ) (

~ since it uses the difference of solutions. By this reason, the problem at hands is underdetermined and, consequently, ill-posed [20,21]. We pose the Inverse Problem (IP) with regularization in order to find a stable and bounded solution. The variational statement [21] with the zero order Tikhonov regularization term is used:

) ( ) )( ( ) ( ) ( , ) ( , ) ( k m jk j m m i k m ik m i j m ij u E u f u A f u A ∆ ∆ + − ∆ − ∆ = α ε .(14)

The first term is a discrepancy of the predictions and observations, the second term poses the zero order Tikhonov regularization, α is the regularization parameter, Ejk is the unite matrix. The regularization term has the form

2 / ) ( 2 / ) ( 2 ) ( 2 ) ( b u u j n j j n j + ∆ = ∆ ∑ ∑ (15)

and ensures the boundedness of b. The minimum of (15) occurs at

m j m n j u u n b ∆ − = ∆ − = ∑ ) ( 1 (16)

So, the minimum attainable error of ∆u (j) is restricted by the mean value:

) ( 1 j m n j m u n u ∆ = ∆ ∑ . (17)

Distance between solutions as the measure of the error

As we mentioned before, the difference between solutions contains some information on errors. Herein, we consider the global (in sense of L2 norm for the grid functions) estimates of errors.

If the relation

2 2 ) 2 ( ) 1 ( 2 ~L h h L h h u u u u − ⋅ ≥ − (18)

holds for numerical solutions u (1) and u (2) , the following contention may be stated.

The norm of approximate solution u (2) error is less than the norm of difference between solutions

2 2 2 , 1 ) 2 ( ~L L du u u ≤ − . (19)

This expression may be easily proved using the triangle inequality [16].

Unfortunately, the information on the errors ordering is not available as a rule.

Fortunately, the additional analysis of distances between solutions may be useful in this situation. For this purpose we should expand the set of analyzed solutions above two. Let u (1) be the maximally inaccurate solution in the ensemble of K numerical solutions. We compare

subsets of distances 2 , 1 L j du and 2 , L j i du ) 1 ( ≠ i . If 2 2 ) ( ) 1 ( ~L i h h L h h u u u u − >> −

(the selected solution is especially inaccurate), the total set of distances between solutions splits into a subset specified by great values of . This situation may be found visually, if the distances between solutions are distributed along a line in accordance with their magnitude. In this situation u (1) may be easily found by the outliers.

The maximum of the distance from zero to maximal value in the cluster 2 , L j i du is assumed to be the upper bound of the cluster δ1 containing distances between "accurate" solutions. The minimum of 2 , 1 L j du is assumed to be a down border of the cluster δ2 containing the distances between "accurate" solutions and most inaccurate one (u (1) ).

The separation of distances between solutions into clusters may be considered as evidence of the existence of solutions with significantly different error norms, that may be stated as the following rule:

If the distance between the clusters is greater than the size of the cluster of accurate solutions

1 1 2 δ δ δ > − , (20) then 2 2 , 1 ) ( ~L i L i du u u ≤ − . (21)

We may use the differences between numerical solutions in different ways since the errors may be of the close magnitudes and the above analysis does not operate. Let's assume these errors ∆u (1) , ∆u (2) to belong hyperspheres with the center at zero point. If the errors are orthogonal, the distance between numerical solutions (hypotenuse) is greater any leg

h k k u u u u u d ) ( ) ( ) 2 ( ) 1 ( 2 , 1 − = ∆ ≥ ∆ − ∆ = (22)

This relation resembles the famous "hypercircle method" [6], unfortunately in certain imprecise form.

The error estimate ( 22) may be naturally extended on the ensemble of K solutions as follows:

2 ) ( max , L k k u d ∆ ≥(23) ) , 1 ( , max )( ) ( max , K i u u d i k i k = − =

The strict orthogonality of approximation errors is not observed in numerical tests [19]. However, the approximation errors are not collinear also. Numerical tests demonstrate the angles between the approximation errors demonstrates α<β and the low boundary may be described as α(β)= β/3. We calculate the angle β using truncation errors δu (j) computed by postprocessor [4] and assume α(β)= β/3 that engenders the estimate:

) 2 , 1 ( , ) 2 / sin( 1 . 1 2 ) ( ) 2 ( ) 1 ( = ∆ > − ⋅ i u u u i α (24)

The comparison of the error estimators

The above considered error estimators are nonintrusive and are based on a postprocessor. We list and discuss the efficiency index (obtained for numerical solution of the compressible Euler equations [14][15][16][17][18][19], containing shock waves and contact discontinuities), the order of the unresolved error and computational expense for this estimators.

Runges' rule

Numerical tests demonstrate the efficiency index Ieff~0.1÷10.

The norm of the unremovable part of the error for this approach has the asymptotics

) ( 1 j h O e = .

This approach uses several consequent grids (minimum two) and, so it is of the medium computational expense.

Richardson extrapolation

Numerical tests [14] demonstrate the Richardson extrapolation to enable the efficiency index Ieff≈1.

The unremovable error is determined by the upper order terms neglected at asymptotic range

) ( 2 j h O e = .

This approach requires four (or greater number) consequent grids and is of the extremely high computational expense.

Inverse Problem

The efficiency index for IP based error estimation is in the range Ieff≈0.25÷4 for K from K=2 and K=13.

The unremovable part of the error e is , where min j is the minimal approximation error over the set of solutions.

This approach requires three (or greater number) independent numerical solutions, obtained on the same grid, and demonstrates the low computational expense.

Distances between solutions

The distance between solutions may be used in several manners that applies the maximum distance between solution (diameter of ensemble), the angle between truncation errors, the analysis of the distances between solutions (the detection of the most imprecise solution).

Diameter of ensemble.

Numerical tests [18] for the openFOAM package demonstrate that the distances between solutions may be used as the error estimators. For the ensemble of five solutions the results of [18] shows Ieff≈0.6÷4. Tests by [19] demonstrates the efficiency index to be in the range Ieff~0.04÷1.5 (K=2) and Ieff~1.1÷1.5, for K=13.

Angle between truncation errors. The account of the angle between the truncation errors (24) enables the estimation that provides the effectivity index in the range Ieff~0.9÷4.52 [19].

The norm of the unremovable part of the error for this approach has the asymptotics

) ( 1 j h O e = .

This approach requires two independent numerical solutions on the same grid and is of the low computational expense.

Analysis of the distances between solutions.

Numerical tests [16] demonstrate the efficiency index for the triangle inequality based estimation (Eqs. ( 19), (21)) is in the range Ieff≈0.75÷2.3.

The norm of the unremovable part of the error for this approach has the asymptotics ) ( 1 j h O e = . This approach requires three (or greater number) independent numerical solutions on the same grid and is of the low computational expense. Unfortunately, it operates only for algorithms having significantly different magnitudes of approximation error. This approach fails for certain sets of solutions ( [18]).

Conclusions

Rigorous a posteriori error estimation for computational Fluid Dynamics at present is practically impossible due to nonlinearity and the discontinuities that may occur and migrate along the flow field. In this situation, several nonstrict (weak) forms of a posteriori estimation of the approximation error may be considered. They either do not provide the error norm estimation in the form of inequalities or provide values of the effectivity index to be less than unit. The best quality of estimates are provided by the Richardson extrapolation, unfortunately for the cost of extremely high computational burden.

However, several new nonstrict forms of a posteriori estimation of the approximation error (based on the ensemble of methods) provide inexpensive estimation of the error norm. The Inverse Problem based error estimation provides the inexpensive estimation of the point-wise error. These approaches hold the greatest promise for the approximation error estimation. These estimators provides the effectivity index