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.
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 [
Au~ = f (1) and corresponding discrete operator
Ahuh = fh (2) that approximates the system on some grid.
In further presentation we consider the numerical solution uh to be a grid function (vector uh∈RM, M is the number of grid points), u~h ∈ R M to be the projection of true solution onto the grid, ∆u~h = uh − u~h 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 u h(i) ∈ R M obtained by independent numerical algorithms (i=1…K) is the number of used algorithm). u h(i) − u~ = ri is the distance between true and h L2 approximate solutions.
Two main options to approximation error exist.
A priori error estimation the estimation
of the
∆u ≤ Chn (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 [
∞ m ∂ m+1u~ δu = ∑ Cm h (4)
m=n
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 [
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 ∂x m+1 ∆u ≤ E(uh ) (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 [
Ie(ffj) = ∆u(i)
h L ∆u~(i) h 2 (6) L2
One may treat the norms of the true error and estimation error as the radii of hyperspheres rexact = ∆u~(i)
h solutions u (i)
h hyperspheres
L2 and rest = ∆uh(i) L . Thus, the numerical
2
with the centre at
are located at surfaces of concentric
u~h and radii ri
(unknown). The relation
I e(fif) ≥ 1
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 [
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 I e(fif) ≥ 1.
We note such error estimators as nonstrict (weak) ones.
The first domain of these approaches forms the defect
correction methods [
Ah (uh )∆uh = δ h . (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 [
Ah (uh )uhrefined = fh +δ h , (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 [
The heuristic rule by C. Runge [
The standards for verification and validation [
The Inverse Problem based approach (IP) [
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 [
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.
From the historic viewpoint the first a posteriori error
indicator is based on the heuristic rule by C. Runge [
The Runge's rule can be easily expanded to the Richardson extrapolation.
The approximation error order that is observed in CFD applications assumes the discrete form
∆u = uh − u~h = C1h j1 + C2h j2 + C3h j3 + .... (9)
where j1, j2, j3, … are positive (sometimes noninteger)
numbers ordered in accordance with the magnitude (for
example, [
The accuracy for the error estimation by Runge's rule has the lowest order O(h j1 ) and remains unresolved.
expansion (9) u (q) = u~ + Chqn
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) [
The pointwise (m is the grid point number) results of numerical computation for three meshes of different steps hq may be presented as:
um(1) = u~m + Ck h1nm um(2) = u~m + Ck h2nm (10)
um(3) = u~m + Ck h3nm .
This system is defined for the most rough grid and
may be resolved regarding u~m , Cm , nm
by several
methods described in [
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
O(h j2 ) and remains unresolved quantitatively that excludes estimates by inequalities of the type (5).
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 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 u~h,m , the approximation error (also unknown) for i-th solution is denoted as ∆um(i) .
The differences of solutions are equal to the differences of the approximation errors and, hence, contain some information regarding the unknown errors ∆um(i) : dij,m = u(i) − u( j) = ∆u(i) − ∆u( j) (11)
m m m m
We treat this information in accordance with the
approach described by [
One may
obtain ∆um(i) .
N = n ⋅ (n − 1) / 2 relations on the set of n numerical solutions:
Aij ∆u( j) = fi,m . (12)
m
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 1 1 0 −1 0 1 0 ∆u m(1) f1,m u m(1) − u m(2) −1 ∆u (2) = f 2,m = u (1) − u (3) . (13)
m m m −1 ∆u m(3) f3,m u m(2) − u m(3)
The solution of these equations is invariant relatively
the transformation ∆u( j) = ∆u~( j) + b since it uses the
m m
difference of solutions. By this reason, the problem at
hands is underdetermined and, consequently, ill-posed
[
j j and ensures the boundedness of b. The minimum of (15) occurs at b = − 1 n
∑ ∆u~m( j) = −∆um (16) n j
So, the minimum attainable error of ∆u(j) is restricted by the mean value: ∆u m = 1 n
∑ ∆u~m( j) . (17) n j
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 u~h − u h(1) ≥ 2 ⋅ u~h − u h(2) (18)
L2 L2 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 u~ − u(2)
L2 2
This expression may be easily proved using the
triangle inequality [
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 ≤ du1,2 L . (19) subsets of distances du1, j L 2 and dui, j L 2 u~h − u h(1) L2 >> u~h − u h(i) L2 (the selected solution is especially inaccurate), the total set of distances between solutions splits into a subset specified by great values of (i ≠ 1) . If du1, j L (distances from accurate solutions to 2 inaccurate one) and a subset of distances between more accurate solutions dui, j L (i ≠ 1) . This situation may 2 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 dui, j L is assumed to be the upper 2 bound of the cluster δ1 containing distances between “accurate” solutions. The minimum of du1, j L 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
δ 2 − δ1 > δ1 , (20) then u~ − u(i)
L2 ≤ du1,i L2 . (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 d1,2 = ∆u (1) − ∆u (2) ≥ ∆u (k ) = u (k ) − u~h (22)
This relation resembles the famous "hypercircle
method" [
The error estimate (22) may be naturally extended on the ensemble of K solutions as follows: d k ,max ≥ ∆u (k )
L2 (23) errors α = arccos dk ,max = max u(k ) − u(i) , (i = 1, K )
i
The strict orthogonality of approximation errors is not
observed in numerical tests [
to be in the range
∆u (1) ⋅ ∆u (2)
30°÷44°. The angles between the truncation errors β are
observed in the range 58°÷64°. Practically all tests
demonstrates α<β and the low boundary may be
described as α(β)= β/3. We calculate the angle β using
truncation errors δu(j) computed by postprocessor [
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 [
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 e = O(h j1 ) .
This approach uses several consequent grids (minimum two) and, so it is of the medium computational expense.
Numerical tests [
The unremovable error is determined by the upper order terms neglected at asymptotic range e = O(h j2 ) .
This approach requires four (or greater number) consequent grids and is of the extremely high computational expense.
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
1 n ∆um = e = ∑ ∆u~m( j)
n j
The norm of the unremovable part of the error for this approach has the asymptotics e = O(h jmin ) , where jmin 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.
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 [
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 [
The norm of the unremovable part of the error for this approach has the asymptotics e = O(h j1 ) .
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 [
The norm of the unremovable part of the error for this approach has the asymptotics e = O(h j1 ) .
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 ([
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 0.3 ≤ Ie(fif) ≤ 5 that may be considered as the acceptable range of the for CFD applications.
This work was supported by grants of RFBR № 1901-00402 and № 20-01-00358.