=Paper=
{{Paper
|id=Vol-2763/CPT2020_paper_s6-3
|storemode=property
|title=Nonstrict methods for a posteriori error estimation
|pdfUrl=https://ceur-ws.org/Vol-2763/CPT2020_paper_s6-3.pdf
|volume=Vol-2763
|authors=A.K. Alekseev,A.E. Bondarev
}}
==Nonstrict methods for a posteriori error estimation==
https://ceur-ws.org/Vol-2763/CPT2020_paper_s6-3.pdf
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.