Discussion: Criticisms of the “Correction Factor” Verification Method 1

Abstract

The technical brief by Roache 1 presents seven criticisms of the correction factor verification method proposed by the authors and colleagues in 234. The authors acknowledge Dr. Roache’s careful review and insightful comments of our work. We have addressed Dr. Roache’s first criticism and revised the correction factor uncertainty estimates; however, we rebut the other criticisms.1. The uncorrected Uk and corrected UkC solution uncertainty estimates given by Eqs. (33) and (34) in 2 are deficient for correction factor Ck⩽1 and Ck=1, respectively, in only providing 50% uncertainty estimate (confidence level), as pointed out by Dr. Roache. In 4, Eq. (33) in 2 was revised for proper behavior for Ck<1: increasing factor of safety for decreasing Ck (i.e., distance from the asymptotic range) similarly as is case for Ck>1. Subsequently, revisions were also made for proper behavior for Ck=1 for both Eqs. (33) and (34) in 2: provision for 10% factor of safety in the limit Ck=1 while smoothly merging with previous correction factor uncertainty estimates for |1−Ck|⩾0.125 (uncorrected solution) and 0.25 (corrected solution) given in 2 and 4. Incorporating both revisions the uncorrected Uk and corrected UkC solution uncertainty estimates are given by: (1)Uk=[9.61−Ck2+1.1]|δREk1*||1−Ck|<0.125[2|1−Ck|+1]|δREk1*||1−Ck|⩾0.125(2)UkC=[2.41−Ck2+0.1]|δREk1*||1−Ck|<0.25[|1−Ck|]|δREk1*||1−Ck|⩾0.25Ck is the correction factor and δREk1* and pk are the one-term Richardson extrapolation (RE) estimates for error and order of accuracy defined by (3)Ck=rkpk−1rkpkest−1(4)δREk1*=εk21rkpk−1(5)pk=lnεk32/εk21lnrkrk is the kth input-parameter (e.g., grid spacing or time step) uniform refinement ratio. εk are solution changes between medium (2) and fine (1), and coarse (3) and medium input parameter values corrected for iterative errors. pkest is an estimate for the limiting order of accuracy of the first term in the RE expansion as spacing size goes to zero and the asymptotic range is reached so that Ck→1. Usually approximation is made that pkest is the theoretical order of the numerical method pkth, e.g., =2 for a second-order accurate method. Replacing equation (5) by (6)pk=lnεk32/εk21lnrk21−1lnrk21 [lnrk32pk−1−lnrk21pk−1]enables use of equations (1) and (2) for nonuniform refinement ratio. Equation (6) corrects equation (31) in 2 for a typographical sign error. 2. The observed order of accuracy is not discarded, but used in defining the correction factor Ck (3), which is then used in defining a corrected one-term RE error estimate =CkδREk1* and uncertainty estimates for uncorrected Uk (1) and corrected Ukc (2) solutions, as shown above and in detail in 2. Correction factors provide a quantitative metric for defining distance of solutions from the asymptotic range and approximately account for the effects of higher-order terms in making error and uncertainty estimates. Correction factors originally based on confirmation studies for 1-D wave 2 and 2-D Laplace equation analytical benchmarks, which showed that the one-term RE error estimate (4) has correct form, but one-term RE order-of-accuracy estimate (5) is poor except in asymptotic range. Multiplication of (4) by Ck (3) provides improved error and uncertainty estimates. For uncorrected solutions, uncertainty estimate Uk (1) based on the absolute value of the corrected error estimate plus the amount of the correction. For corrected solutions (i.e., corrected error estimate is used both in sign and magnitude to define numerical benchmark Sc=S−CkδREk1*),UkC (2) based on the absolute value of the amount of the correction. 3. Reference 4 shows that the correction factor approach is equivalent to the GCI, but with a variable factor of safety (FS), which increases with distance of solutions from the asymptotic range. Ck (3) provides metric for estimating distance of solutions from the asymptotic range: =1 when solutions are in asymptotic range; <1 when pk<pkest; and >1 when pk>pkest. For GCI approach, FS is constant for all Ck: 1.25 for careful grid studies otherwise 3. For Ck approach, FS varies linearly with Ck with slope 2 (uncorrected solution) and 1 (corrected solution) and symmetric about Ck=1. The intersection points between Ck and GCI approaches depends on value FS used in GCI, e.g., for FS=1.25 intersection points are Ck=0.875,1.125 and (0.75,1.25) for uncorrected and corrected solutions, respectively. When solutions are between the intersections points (closer to the asymptotic range), GCI approach is more conservative than Ck approach. When solutions are outside the intersection points (further from the asymptotic range), GCI approach is less conservative than Ck approach. The previously mentioned analytical benchmarks confirmed the FS slope predicted by correction factors and admittedly may not be best for all cases (an interesting topic for future research). Variable FS that increases with distance from the asymptotic range is a “common-sense” advantage of correction factor approach compared to GCI, especially for practical applications where solutions are often further from the asymptotic range. In retrospect, correction factor approach has similarities to that proposed by Celik and Karatekin 5, who in present notation propose UC&K=|CkδRE*| (equivalent equation (9a) of 5). UC&K is similar to correction factors for Ck>1 in providing variable FS which increases with distance from the asymptotic range; however, less conservative than correction factors which includes the amount of the correction. Also, for Ck=0 UC&K suffers 50% uncertainty and for Ck<1 gives an unacceptable result of FS<1. Figure 1 compares FS predicted by correction factors, GCI, and UC&K.4. Our claim of GCI being more conservative for solutions closer to asymptotic range and more importantly less conservative for solutions further from the asymptotic range in comparison to correction factors is based on our comparison of factor of safety for both approaches, as discussed above as well as comparisons between two approaches for analytical and numerical benchmarks and practical applications. The cases in 34 are practical applications with solutions further from the asymptotic range, especially 3; therefore, correction factor uncertainties are more conservative than GCI. For the Series 60 model ship 3, verification results were obtained on four grids wherein for the finer three grids pG=4.4 CG=3.7 and for the coarser three grids pG=2.3 CG=1.3 for the total resistance coefficient. In the former case, the GCI approach with FS=1.25 estimates uncertainties which are smaller than the correction factor approach by a factor of roughly five. If FS=3 is used, the GCI approach estimates uncertainties which are smaller than the correction factor approach by a factor of roughly 2. For the 5415 surface combatant 4, verification results were obtained on three grids and both correction factor and GCI give similar results for the total resistance coefficient since Ck=0.8 is near the intersection of the two approaches. Ebert and Gorski 6 showed that the correction factor approach is reasonable for pk>pkest, but not for pk<pkest where Ck uncertainty estimates indicated insufficient conservatism. Present revisions no longer have this deficiency, as previously discussed. For analytical benchmarks, as already mentioned, 1-D wave and 2-D Laplace equations originally used to test verification methods. Analytical benchmarks are ideal for confirmation of verification procedures since the exact solution known analytically and modeling errors are zero such that numerical errors can be determined exactly; however, analytical benchmarks are restricted to simple problems. For the 1-D wave equation, verification results obtained on 10 grids with Ck values for maximum wave amplitude in the range 0.5⩽Ck⩽1. For the uncorrected simulation, both Uk and GCI (with FS=1.25) banded the true error with the latter more conservative. The corrected simulation, UkC banded the true error, whereas GCI (extended for use with corrected solutions) failed to band the true error for the two coarsest grids. Verification results on 6 grids for the 2-D Laplace equation with constant Dirichlet boundary conditions showed Ck values at a fixed point in the range 0.8⩽Ck⩽1.0. For the uncorrected simulations, both Uk and GCI (with FS=1.25) banded the true error with the latter more conservative for solutions closer to the asymptotic range and less conservative for solutions further from the asymptotic range. For the corrected simulation, both UkC and GCI banded the true error. More recently, verification results were obtained for the Blasius boundary layer analytical benchmark 7, on 5 grids with Ck values for the axial velocity profile at quarter-plate in the range 0.75<Ck<5.0. For the finest grid triplet, both Ck and GCI approaches band the true error, but for the coarser grid triplets, the GCI approach failed to band the true error for the inner part of the boundary layer. Results also showed difficulties for complex nonlinear analytical benchmarks such as the necessity of restricting the analysis to the region of flow for which boundary layer theory is valid. In preparation of the present discussion, additional comparison of verification procedures used analytical benchmark data for 2-D Laplace equation with linear/logarithmic varying Dirichlet boundary conditions from 8 and numerical benchmark data for k-ε model turbulent flow backward facing step from 5. Numerical benchmarks use grid-independent solutions to estimate the numerical benchmark SC and numerical error δSN* as the difference between the fine Sk1 and grid independent solutions δSN*=Sk1−Sc. Numerical benchmarks are not restricted to simple problems; however, are less-ideal than analytical benchmarks for confirmation of verification procedures since assumptions must be made that benchmark solution is without numerical error and modeling errors are same for all solutions both of which are suspect for complex problems. Results for the former are similar to previously discussed 2-D Laplace equation results. Note that for both 1-D wave and 2-D Laplace equation cases the UC&K approach fails to bound the true error since Ck<1 for reasons discussed previously. Results for the latter, show that for the uncorrected solution both correction factor and UC&K approaches with pkest=1 bound the estimated error with correction factors more conservative, whereas the GCI approach (with FS=1.25) fails to bound the estimated error (based on previously published “grid independent” solutions with same turbulence model) for the finest and a medium grid. For the corrected simulation, both correction factor and GCI approaches fail to band the estimated error. Results were also obtained for the correction factor approach with pkest=2; since, 5 uses the hybrid numerical method with variable pth between 1 and 2. For the uncorrected simulation, correction factor trends are similar to GCI, but more conservative. For the corrected simulation, the correction factor approach bands the estimated error, except on the finest grid. In 5, UC&K and GCI (with FS=3) are compared showing latter more conservative than former; however, for such careful grid study GCI with FS=1.25 is more appropriate. Data and figures showing comparisons of correction factor, GCI, and Celik and Karatekin verification methods for aforementioned analytical and numerical benchmarks are available from the authors upon request. The overall results show that when solutions are closer to the asymptotic range both correction factor and GCI approaches provide reasonable results, although the GCI approach may be over conservative in comparison to the correction factor approach. When solutions are further from the asymptotic range, as is often the case with practical applications, correction factors have the advantage of increased factor of safety and GCI approach may be under conservative in comparison to the correction factor approach. The Celik and Karatekin 5 approach is reasonable for Ck>1 although less conservative than correction factor approach and not reasonable for Ck⩽1.5. Strictly speaking, the corrected solution (i.e., simulation corrected using deterministic estimate of sign and magnitude of numerical error) will not satisfy the same conservation properties (e.g., mass and momentum) as the uncorrected solution. However, if solutions close to the asymptotic range one would expect the correction and the lack of conservation for the numerical benchmark solution to be small. Situations that might prevent correction of the solution include variability in the observed order of accuracy, lack of complete iterative convergence, and solutions further from the asymptotic range. Situations also exist when it is useful to make use of the corrected solution. Confirmation of verification procedures using numerical benchmarks in which the corrected solution is assumed without numerical error and modeling errors are assumed same for all solutions 9 and for practical applications as aid in assessing modeling errors when monotonic convergence established for multiple grid triplets and/or iterative convergence and resource issues for very fine grids 34. Furthermore, the concept of deterministic error estimate for simulations seems appropriate and been advocated by others. We agree with the caveat that it is only useful as a deterministic estimate under certain circumstances (i.e., solutions sufficiently close to the asymptotic range) and in this case, the uncertainty estimate based on an estimate of the error in that estimate. We never claimed using the corrected solution addresses the criticism that the numerical error is deterministic and not stochastic. We disagree with Dr. Roache’s formulation for estimating the uncertainty of corrected solutions; since, the amount of conservatism is a function of grid refinement ratio and for rk<2 gives the unacceptable result that the uncertainty in the corrected solution is larger than the error in the uncorrected solution. Therefore, we recommend the revised correction factor uncertainty estimates for the corrected solution UkC (2). The revised uncertainty estimates given by equation (2) may improve the Eca and Hoekstra 8 results showing that uncertainty estimates for the corrected solution from the correction factor approach do not bound the error for the 2-D Laplace equation analytical benchmark. The data presented in 8 are insufficient to reproduce all their results. However, we are able to reproduce a subset of their results, which show that correction factor corrected solution uncertainty estimates bound the error, as we have also shown for the 2-D Laplace equation. 6. We agree that analytical benchmarks are useful to confirm verification procedures; however, we disagree that the ensemble of problems in 10 or data of 9 provide statistical evidence for establishing 95% confidence level. We find no statistical distributions in 9 or 10 while 10 assumes normal distributions with the claim that twice the standard deviation leads to a 95% confidence level with no supporting evidence. Truncation errors are systematic errors with strong spatial and temporal correlations; therefore, we do not expect errors for single problems (i.e., individual user, code, model, grid-type, etc.) to display normal distributions. One must consider issues of replication level, as with experimental uncertainty analysis. Recently, authors have proposed a method for establishing probabilistic confidence intervals for CFD codes 11 using N-version testing 12 wherein multiple codes or users, models, grid types etc. for specific benchmark applications are used to establish normal distributions. 7. We do not understand how you can accuse us of ignoring issues related to erroneous identification of monotonic convergence when we specifically provided discussion in 234, including reference to 13. We have also extensively discussed in 234 a minimum of three solutions are required for verification studies and desirability and issues for obtaining more than three solutions, including reference to 8. We are indebted to Dr. Roache for pointing out the possibility of the oscillatory divergence convergence condition and to Professor Celik for pointing out his discussion in 5 of oscillatory convergence.