Thursday, September 22, 2011

Mystery surrounding order of accuracy in CFD...

Just came back from a trip to Beijing, where I gave two lectures in a Workshop, one on high-order CFD methods, and the other on the quadrature method of moment for multiphase flows. On the last day, I got to see the beautiful Great Wall at Mu Tian Yu. See the top picture for a magnificent view! Personally I believe this is a better tourist attraction than the more famous Great Wall at Ba Da Lin because it is less crowded. I am now paying a price for the trip: I am still recovering from the jet lag.

In the Workshop at Beijing as well as in many other technical conferences, I often heard people talking about and making mistakes on order of accuracy (OOA) in CFD. As often as it is mentioned, there appears to be as much confusion about it too. There are two major mistakes regarding OOA:
  1. If a flow solver produces a more accurate solution than another, it must posses  higher order accuracy. 
  2. The OOA is the same as the order of the polynomial approximating the (unknown) solution.

I still remember an incident in an AIAA conference (maybe at Reno). A guy believed like a religion that the OOA is the same as the order of the reconstructed solution polynomial in a finite volume scheme. In fact, he accused the finite volume community of lying about the accuracy of their favorite schemes. You cannot help feeling bad for the guy with such a strong (but wrong) conviction. The truth is that the OOA of a finite volume scheme with a degree p solution polynomial is (p+1).

Let's define OOA. A scheme is pth order if the solution error is proportional to (dx)^p, where dx is the mesh size, i.e.,

Error = C(dx)^p,

where C is a constant depending on the particular scheme. Or equivalently,

log(Error) = p*log(dx) + log(C).

In other words, p is the slope of the dx-Error plot in a log-log scale. Strictly speaking, OOA has nothing to do with the magnitude of the solution error, but only related to the rate of error reduction with mesh refinement, i.e., the slope.

How do you actually measure the OOA? If you know the exact solution, you need only perform two simulations: one on a coarse mesh, and the other on a fine mesh with half the mesh size of the coarse mesh. With a simple derivation, you can obtain the order p from

Error_(dx)/Error_(dx/2) = 2^p, or p = log_2[Error_(dx)/Error_(dx/2)]

where Error_(dx) is the error on the coarse mesh, and  Error_(dx/2) the error on the fine mesh.

More later on high-order methods...


  1. Hello Dr. Wang.

    Does not the OOA of a finite volume scheme depends not only on the order of the reconstructed solution polynomial but also on the order of the integration method used to carry out the integral over the control volume faces? Most finite volume schemes assume that the value of what is being integrated over the control volume is a constant, therefore such schemes will have an OOA no greater than two, irrespective of the order of the reconstructed solution polynomial.


  2. Hi Frank,

    Thanks for the comments. I completely agree with you that the precision of the surface flux integration is as important as the solution polynomial in determining the OOA of the finite volume method. Normally, Gauss quadratures are used for such integrations, which must be exact for a degree p polynomial in order to achieve a (p+1)th order accuracy. The mid-point integration rule is fine for both piece-wise constant and linear solution polynomials.


  3. Dear Dr. Wang,
    I'a a PhD candidate, writing the thesis on FV simulations of spinning detonations using my in-house reactive Euler code. I've implemented AMR methods on unstructured grids, and now I am facing the problem of assessing the accuracy and efficiency of various refinement citeria. I am a little confused, since I am not sure if the method of determining the order of accuracy, as the one presented by you above, can be applied to unsteady problems.
    I've started with Sod's test (from book of Toro), testing up to 10 levels of fully refined, regular grid and I've found for my first order HLLC solver, that the p(slope) tends to lower as the grid is being refined. (form 0.5 for 1st level to 0.1 for 10-th level).
    When I use 10-th level solution as a reference solution for L1 error norm, instead of accurate solution, the slope rises form 1 to 2.7 as the grid is being refined.
    Any of the above results doesn't make any sense to me, so I've started to look for explanation. So my question is
    Is the method of determining the order of accuracy presented above, valid for unsteady solutions of nonlinear problems? I believe, the nonlinearity of the problem may cause the initial error to develop in time, causing for example wrong wave speed estimates and finally displacement of discontinuity in Sod's problem.
    Therefore my question is, maybe such comparisons should be done not for the whole simulations, but for single "coarse" time step, like in Richardson's method?

    Thank you in advance for reply,
    Michal Folusiak
    Warsaw University of Technology, Poland

  4. Michal,

    I suggest that you study the error at a certain time, say, time=1 second, to assess the order of accuracy for unsteady problems. Obviously, the error at that time depends on both spatial and time discretization methods. You could assess the spatial error by using very small time steps.

    Strictly speaking, you can only obtain the designed order of accuracy if your solution is sufficiently smooth. I therefore recommend that you use a smooth analytical solution to assess your Euler solver first as a verification case. Then you can extend the same approach to study the numerical OOA for no smooth

  5. for non smooth solution. However, if you failed to obtain the designed order of accuracy, you can always blame the non smooth solution!

    Good luck,