Nonstrict methods for a posteriori error estimation A.K. Alekseev1, A.E. Bondarev1 aleksey.k.alekseev@gmail.com | bond@keldysh.ru Keldysh Institute of Applied Mathematics RAS, Moscow, Russia; The paper is devoted to comparison of a posteriori methods (based on the precomputed solutions) for approximation error estimation. 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. We pay the special attention to the nonstrict methods that either cannot be presented in a form of inequalities, or demonstrate the effectivity index of an estimator to be below unit. Several new, computationally inexpensive methods for both the point-wise error and the error norm estimation are considered. They are nonintrusive, realized by postprocessing and provide a successful compromise of the reliability and computational efforts. Methods based on the use of an ensemble of independent solutions can be implemented by constructing a generalized computational experiment, which sharply increases the speed and efficiency of the assessment. Keywords: approximation error, partial differential equations, a posteriori error estimation, Richardson extrapolation, Runge rule, nonstrict estimators. (source term of differential approximation [3]) that may 1. Introduction be formally expressed as The approximation error is omnipresent at the ∞ ∂ m +1u~ numerical solutions of Partial Differential Equations δu = ∑ C m h m (4) (PDE) due to the discretization at numerical statements. m=n ∂x m +1 The error estimation is of the high current interest in view for PDE systems of the first order. This expression of the need for the verification of software and numerical contains the infinite number of terms so, the first (minor solutions. For example, the corresponding issues are order) term is commonly used. It may be computed by stated in Computational Fluid Dynamics (CFD) in the many ways including the special postprocessor [4]. form of standards [1,2]. Let's discuss main approaches to A priori error estimates have an universal form. the estimation of the approximation error. We consider a Unfortunately, the unknown constant prevents it from to PDE system in the operator form be used in applications. Au~ = f (1) A posteriori error estimator usually has the form and corresponding discrete operator ∆u ≤ E (uh ) (5) Ah uh = f h (2) and is determined by some computable error indicator E(uh). This estimator depends on the previously that approximates the system on some grid. computed numerical solution uh and, thus, has a minor In further presentation we consider the numerical generality. Fortunately, it may be applied to practical solution uh to be a grid function (vector uh∈RM, M is the computations since has no unknown constants. ~ ∈R number of grid points), u M to be the projection of The highly efficient technique is developed for a h ~ = u − u~ to be the true true solution onto the grid, ∆u posteriori error estimation in the domain of the finite- h h h element analysis [5,7,8]. In accordance with [5], the approximation error, ∆uh to be some estimate of this quality of a posteriori error estimation may be expressed error. L2 - based norm is used for a comparison of these via the effectivity index of estimator that is equal to the vectors. We may also use the set of numerical solutions relation of the estimated error norm to the true error u h(i ) ∈ R M obtained by independent numerical norm: algorithms (i=1…K) is the number of used algorithm). ∆uh(i ) u h(i ) − u~h = ri is the distance between true and I ( j) eff = ~ (i ) L2 (6) L2 ∆u h L2 approximate solutions. One may treat the norms of the true error and Two main options to the estimation of the estimation error as the radii of hyperspheres approximation error exist. A priori error estimation rexact = ∆u~h(i ) and rest = ∆uh (i ) . Thus, the numerical L2 L2 n ∆u ≤ Ch (3) solutions u h(i ) are located at surfaces of concentric is commonly used at the design and the theoretical ~ and radii ri analysis for the determination of the convergence order. hyperspheres with the centre at u h Herein h is the step of discretization, n is the order of (unknown). The relation I (i ) ≥ 1 means that the eff approximation, C is an unknown constant. This estimation is usually related to the truncation error hypersphere, containing the true error, belongs to the hypersphere defined by the estimator. So, in order to Copyright © 2020 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0) provide the reliable estimation, the effectivity index Thus, the computationally inexpensive nonstrict a should be greater the unit. On the other hand, the posteriori estimation of approximation error is of the estimation should be not too pessimistic, so the value of major interest in CFD from the viewpoint of verification the effectivity index should be not too great. For the of codes and solutions. The simultaneous use of several finite element methods, used in the domain of elliptic nonstrict methods may have some prospects from the equations (usually engendering highly regular solutions), viewpoint of reliability increasing. (i ) the acceptable range, according [5], is 1 ≤ I eff ≤ 3 . 2. Runge rule However, the boundaries of this inequality seem to be dependent on the problem at hand. Numerical tests for From the historic viewpoint the first a posteriori error CFD domain demonstrate the efficiency index to belong indicator is based on the heuristic rule by C. Runge [5]. If the range (0.3,5). For the nonlinear problems containing the difference between two approximate solutions discontinuities (that is common case for CFD) the computed on a coarse mesh uh and the refined mesh uh,ref progress of a posteriori error estimation in the rigorous is small, then both are assumed to be close to the exact solution. The Runge’s rule can be considered as the first form of inequality (5) is limited. As an alternative, some less rigourous methods are a posteriori error indicator |ε(uh)- ε(uh,ref)|=ERunge(uh- uh,ref) if one uses certain functional of the flow variables. employed. These methods provide the estimation of ∆uh It is the basis for the stopping criterion by merging of without any strict inequality. Significant number of (i ) some functional at the mesh refinement. However, such estimators do not met the condition I eff ≥ 1 . relations do not guarantee convergence of the total We note such error estimators as nonstrict (weak) solution or other valuable functionals. From a practical ones. needs perspective one should desire the quantitative The first domain of these approaches forms the defect estimate of the form uh − u~ ≤ δ with computable δ . correction methods [9-11]. Some part of these methods The Runge's rule can be easily expanded to the [9,10] apply some approximation of the truncation error Richardson extrapolation. δuh in order to disturb the main system. The additional The approximation error order that is observed in equation for the error transformation occurs CFD applications assumes the discrete form Ah (u h )∆u h = δ h . (7) ∆u = u h − u~h = C1h j1 + C 2 h j2 + C3 h j3 + .... (9) These methods are rather laborious since imply where j1, j2, j3, … are positive (sometimes noninteger) coding, debugging and solving of an additional problem. numbers ordered in accordance with the magnitude (for A bit less laborious version of defect correction example, [12]). methods [11] have the appearance The accuracy for the error estimation by Runge's rule Ah (u h )u hrefined = f h + δ h , (8) j has the lowest order O (h 1 ) and remains unresolved. that imply the disturbing of the main problem by the source term, which approximates the truncation error. 3. Richardson extrapolation Another branch of nonstrict a posteriori error estimation methods has a non-intrusive form of certain Richardson extrapolation is based on the first term of postprocessor that significantly reduces efforts for coding expansion (9) u ( q ) = u~ + Chqn and operates if the and debugging. It may be based on the Runge rule [5], asymptotic range is achieved for several grids hq (Ck, n Richardson extrapolation (RE) [13,14], Inverse Problem are assumed to be constant that should be verified based approach (IP) [15] or ensemble based methods numerically by expanding the set of grids) [13,14]. For (EM) [16-19]. CFD problems containing shock waves and contact lines The heuristic rule by C. Runge [5] is the basis of [12] the error order is not constant over the flowfield and commonly used stopping criterion by merging of depends on the type of flow structure. So, one should to solutions at the mesh refinement. extend RE for estimation of the local order of The standards for verification and validation [1,2] convergence. recommend the Richardson extrapolation (RE) as the The pointwise (m is the grid point number) results of main tool for the verification. Richardson extrapolation numerical computation for three meshes of different steps provides the pointwise approximation of the error field, hq may be presented as: unfortunately, at the cost of the high computational burden [13,14]. RE provides some generalization of the um(1) = u~m + Ck h1n m Runge's rule. u ( 2 ) = u~ + C h n m (10) m m k 2 The Inverse Problem based approach (IP) [15] enables the pointwise information on the error. u ( 3) m = u~m + Ck h3n m . The computationally cheap approach to a posteriori This system is defined for the most rough grid and error estimation that is based on the ensemble of numerical solutions obtained by independent methods is may be resolved regarding u~m , C m , nm by several offered by [16-19]. However, these approaches do not methods described in [13, 14]. The relations (10) are provide the pointwise information on the approximation valid, if Cm is independent on h and higher order terms error. 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 parameter, Ejk is the unite matrix. The regularization term these grids, the additional refinement is necessary. has the form So, the Richardson extrapolation provides the n n pointwise approximation for the error field at the cost of ∑ (∆u ( j ) )2 / 2 = ∑ (∆u~ ( j ) + b)2 / 2 (15) j j the high computational burden. The accuracy of RE for the error estimation has the and ensures the boundedness of b. The minimum of (15) j occurs at appearance O(h 2 ) and remains unresolved 1 n quantitatively that excludes estimates by inequalities of b=− ∑ ∆u~m( j ) = −∆um (16) the type (5). n j So, the minimum attainable error of ∆u(j) is restricted 4. Approximation error estimation using by the mean value: Inverse Problems 1 n The nonstrict (weak) form of the approximation error ∆u m = ∑ ∆u~m( j ) . (17) n j 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 5. Distance between solutions as the measure consider an ensemble of numerical solutions of the error um(i ) = u~h , m + ∆um(i ) , obtained by K independent As we mentioned before, the difference between solutions contains some information on errors. Herein, algorithms (of different structure, for example, of we consider the global (in sense of L2 norm for the grid different approximation order) on the same grid. The functions) estimates of errors. projection of an exact (unknown) solution on the grid ~ , the approximation error (also If the relation point m is denoted as u h, m u~h − u h(1) ≥ 2 ⋅ u~h − u h( 2 ) (18) (i ) L2 L2 unknown) for i-th solution is denoted as ∆u m . holds for numerical solutions u(1) and u(2), the following The differences of solutions are equal to the contention may be stated. differences of the approximation errors and, hence, The norm of approximate solution u(2) error is less contain some information regarding the unknown errors than the norm of difference between solutions ∆um(i ) : u~ − u ( 2 ) ≤ du1, 2 L . (19) L2 d ij , m = um(i ) − um( j ) = ∆um(i ) − ∆um( j ) (11) 2 This expression may be easily proved using the triangle inequality [16]. We treat this information in accordance with the Unfortunately, the information on the errors ordering approach described by [15] in order to determine the is not available as a rule. approximation error ∆um(i ) . One may obtain Fortunately, the additional analysis of distances between solutions may be useful in this situation. For this N = n ⋅ (n − 1) / 2 relations on the set of n numerical purpose we should expand the set of analyzed solutions solutions: above two. Let u(1) be the maximally inaccurate solution Aij ∆um( j ) = f i , m . (12) in the ensemble of K numerical solutions. We compare The summation over the repeating indexes is used subsets of distances du1, j L2 and dui , j L2 (i ≠ 1) . If elsewhere from this point. For the minimum set of data (three solution) equation (12) assumes the form u~h − u h(1) >> u~h − u h(i ) (the selected solution is L2 L2  1 − 1 0  ∆u   f1, m   u − u  (1) (1) ( 2) especially inaccurate), the total set of distances between   m     m m  ( 2) (1) ( 3) solutions splits into a subset specified by great values of  1 0 − 1 ∆u m  =  f 2, m  =  u − u m m  . (13)  0 1 − 1 ∆u ( 3)   f  u − u ( 2) ( 3)  du1, j (distances from accurate solutions to   m   3, m   m m  L2 The solution of these equations is invariant relatively inaccurate one) and a subset of distances between more the transformation ∆u ( j) m = ∆u~m( j ) + b since it uses the accurate solutions dui , j L2 (i ≠ 1) . This situation may difference of solutions. By this reason, the problem at be found visually, if the distances between solutions are hands is underdetermined and, consequently, ill-posed distributed along a line in accordance with their [20,21]. We pose the Inverse Problem (IP) with magnitude. In this situation u(1) may be easily found by regularization in order to find a stable and bounded the outliers. solution. The variational statement [21] with the zero The maximum of the distance from zero to maximal order Tikhonov regularization term is used: value in the cluster dui , j is assumed to be the upper ε = ( Aij ∆um( j ) − f i ,m )( Aik ∆um( k ) − f i ,m ) + α (∆um( j ) E jk ∆um( k ) ) . (14) L2 The first term is a discrepancy of the predictions and bound of the cluster δ1 containing distances between observations, the second term poses the zero order “accurate” solutions. The minimum of du1, j is Tikhonov regularization, α is the regularization L2 assumed to be a down border of the cluster δ2 containing The norm of the unremovable part of the error for this the distances between “accurate” solutions and most approach has the asymptotics e = O(h j1 ) . inaccurate one (u(1)). The separation of distances between solutions into This approach uses several consequent grids clusters may be considered as evidence of the existence (minimum two) and, so it is of the medium computational of solutions with significantly different error norms, that expense. may be stated as the following rule: Richardson extrapolation If the distance between the clusters is greater than the size of the cluster of accurate solutions Numerical tests [14] demonstrate the Richardson δ 2 − δ1 > δ1 , (20) extrapolation to enable the efficiency index Ieff≈1. The unremovable error is determined by the upper then order terms neglected at asymptotic range u~ − u (i ) ≤ du1,i L . (21) e = O(h j2 ) . L2 2 We may use the differences between numerical This approach requires four (or greater number) solutions in different ways since the errors may be of the consequent grids and is of the extremely high close magnitudes and the above analysis does not computational expense. operate. Let's assume these errors ∆u(1), ∆u(2) to belong hyperspheres with the center at zero point. If the errors Inverse Problem are orthogonal, the distance between numerical solutions The efficiency index for IP based error estimation is (hypotenuse) is greater any leg in the range Ieff≈0.25÷4 for K from K=2 and K=13. d1, 2 = ∆u (1) − ∆u ( 2 ) ≥ ∆u ( k ) = u ( k ) − u~h (22) The unremovable part of the error e is 1 n This relation resembles the famous "hypercircle ∆u m = e = ∑ ∆u~m( j ) method" [6], unfortunately in certain imprecise form. n j The error estimate (22) may be naturally extended on the ensemble of K solutions as follows: The norm of the unremovable part of the error for this d k , max ≥ ∆u ( k ) (23) approach has the asymptotics e = O(h jmin ) , where L2 jmin is the minimal approximation error over the set of (k ) (i ) d k ,max = max u −u , (i = 1, K ) solutions. i This approach requires three (or greater number) The strict orthogonality of approximation errors is not independent numerical solutions, obtained on the same observed in numerical tests [19]. However, the grid, and demonstrates the low computational expense. approximation errors are not collinear also. Numerical tests demonstrate the angles between the approximation Distances between solutions (1) ( 2) errors α = arccos (∆u , ∆u ) to be in the range The distance between solutions may be used in ∆u (1) ⋅ ∆u ( 2 ) several manners that applies the maximum distance 30°÷44°. The angles between the truncation errors β are between solution (diameter of ensemble), the angle observed in the range 58°÷64°. Practically all tests between truncation errors, the analysis of the distances demonstrates α<β and the low boundary may be between solutions (the detection of the most imprecise described as α(β)= β/3. We calculate the angle β using solution). truncation errors δu(j) computed by postprocessor [4] and Diameter of ensemble. assume α(β)= β/3 that engenders the estimate: Numerical tests [18] for the openFOAM package demonstrate that the distances between solutions may be u (1) − u ( 2 ) used as the error estimators. For the ensemble of five 1.1 ⋅ > ∆u (i ) , (i = 1,2) (24) solutions the results of [18] shows Ieff≈0.6÷4. Tests by sin(α / 2) 2 [19] demonstrates the efficiency index to be in the range 6. The comparison of the error estimators Ieff~0.04÷1.5 (K=2) and Ieff~1.1÷1.5, for K=13. Angle between truncation errors. The above considered error estimators are The account of the angle between the truncation nonintrusive and are based on a postprocessor. We list errors (24) enables the estimation that provides the and discuss the efficiency index (obtained for numerical effectivity index in the range Ieff~0.9÷4.52 [19]. solution of the compressible Euler equations [14-19], The norm of the unremovable part of the error for this containing shock waves and contact discontinuities), the approach has the asymptotics e = O(h j1 ) . order of the unresolved error and computational expense for this estimators. This approach requires two independent numerical solutions on the same grid and is of the low Runges' rule computational expense. Analysis of the distances between solutions. Numerical tests demonstrate the efficiency index Ieff~0.1÷10. Numerical tests [16] demonstrate the efficiency index [7] I. Babuska and W. Rheinboldt. A posteriori error for the triangle inequality based estimation (Eqs. (19), estimates for the finite element method. Int. J. Numer. (21)) is in the range Ieff≈0.75÷2.3. Methods Eng. 12: 1597–1615 The norm of the unremovable part of the error for this [8] M Ainsworth. and J. T. Oden, A Posteriori Error Estimation in Finite Element Analysis. Wiley – approach has the asymptotics e = O(h j1 ) . Interscience, NY. (2000). This approach requires three (or greater number) [9] T. Linss and N. Kopteva, A Posteriori Error Estimation independent numerical solutions on the same grid and is for a Defect-Correction Method Applied to Convection- of the low computational expense. Unfortunately, it Diffusion Problems, Int. J. of Numerical Analysis and operates only for algorithms having significantly Modeling, V. 1, N. 1, (2009) 1–16. different magnitudes of approximation error. This [10] J. W. Banks, J. A. F. Hittinger, C. S. Woodward, approach fails for certain sets of solutions ([18]). Numerical error estimation for nonlinear hyperbolic PDEs via nonlinear error transport, CMAME, 213 7. Conclusions (2012) 1-15. [11] Christopher J. Roy, and Anil Raju, Estimation of Rigorous a posteriori error estimation for Discretization Errors Using the Method of Nearby computational Fluid Dynamics at present is practically Problems, AIAA JOURNAL, Vol. 45, No. 6, June 2007 impossible due to nonlinearity and the discontinuities that p 1232-1243 may occur and migrate along the flow field. In this [12] J. W. Banks, T. D. Aslam, Richardson Extrapolation for situation, several nonstrict (weak) forms of a posteriori Linearly Degenerate Discontinuities, Journal of estimation of the approximation error may be considered. Scientific Computing, May 24, 2012 P. 1-15 They either do not provide the error norm estimation in [13] Ch. J. Roy, Grid Convergence Error Analysis for Mixed the form of inequalities or provide values of the –Order Numerical Schemes, AIAA Journal, V. 41, N. 4, effectivity index to be less than unit. The best quality of (2003) 595-604. estimates are provided by the Richardson extrapolation, [14] Alexeev, A.K., Bondarev, A.E.: On Some Features of unfortunately for the cost of extremely high Richardson Extrapolation for Compressible Inviscid computational burden. Flows. Mathematica Montisnigri XL, 42–54 (2017). [15] Alekseev A.K., Bondarev A. E., Kuvshinnikov A. E., A However, several new nonstrict forms of a posteriori posteriori error estimation via differences of numerical estimation of the approximation error (based on the solutions, ICCS 2020. ensemble of methods) provide inexpensive estimation of [16] Alekseev, A.K., Bondarev, A.E., Navon, I.M.: On the error norm. The Inverse Problem based error Triangle Inequality Based Approximation Error estimation provides the inexpensive estimation of the Estimation. arXiv:1708.04604 [physics.comp-ph], point-wise error. These approaches hold the greatest August 16, 2017. promise for the approximation error estimation. These [17] Alekseev A.K., Bondarev A. E., Kuvshinnikov A. E.: (i ) estimators provides the effectivity index 0.3 ≤ I eff ≤ 5 Verification on the Ensemble of Independent Numerical Solutions, In: Rodrigues J. et al. (eds) Computational that may be considered as the acceptable range of the for Science – ICCS 2019. ICCS 2019. Lecture Notes in CFD applications. Computer Science, Springer, Cham, 11540, 315–324 (2019). Acknowledgments [18] Alekseev A.K., Bondarev A. E., Kuvshinnikov A. E., This work was supported by grants of RFBR № 19- On uncertainty quantification via the ensemble of 01-00402 and № 20-01-00358. independent numerical solutions // Journal of Computational Science 42 (2020) 101114, DOI: References 10.1016/j.jocs.2020.101114 [19] A. K. Alexeev, A. E. Bondarev, The features of the of [1] Guide for the Verification and Validation of truncation and approximation errors’ geometry on the Computational Fluid Dynamics Simulations, American ensemble of numerical solutions, KIAM Preprint, Institute of Aeronautics and Astronautics, AIAA-G- Moscow (2019) № 107 (in Russian), DOI: 077-1998, Reston, VA, 1998. 10.20948/prepr-2019-107 [2] Standard for Verification and Validation in [20] Tikhonov, A.N., Arsenin, V.Y.: Solutions of Ill-Posed Computational Fluid Dynamics and Heat Transfer, Problems. Winston and Sons, Washington DC (1977). ASME V&V 20-2009, 2009 [21] Alifanov, O.M., Artyukhin, E.A., Rumyantsev S.V.: [3] Yu.I Shokin, Method of differential approximation. Extreme Methods for Solving Ill-Posed Problems with Springer-Verlag, (1983). Applications to Inverse Heat Transfer Problems. Begell [4] A.K. Alekseev, I.M. Navon, A Posteriori Error House (1995). Estimation by Postprocessor Independent of Flowfield Calculation Method, Computers & Mathematics with About the authors Applications, v. 51, (2006) 397-404. [5] Repin, S.I.: A posteriori estimates for partial differential Alekseev Aleksey K., senior researcher, Doctor of Physical equations. Vol. 4. Walter de Gruyter (2008). and Mathematical Sciences, Keldysh Institute of Applied [6] W. Prager and J. L. Synge. Approximation in elasticity Mathematics RAS, Moscow. E-mail: based on the concept of function spaces, Quart. Appl. aleksey.k.alekseev@gmail.com Math. 5 (1947) 241-269 Bondarev Alexander E., PhD, senior researcher, Keldysh Institute of Applied Mathematics RAS. E-mail: bond@keldysh.ru.