how to determine and prove the order of our cfd code?
 

Theoretically it's possible to prove that our discretization schemes are 2nd order. But during coding, it's possible to introduce errors that make it 1st order.


If possible, I would like to prove for my codes, that they're indeed as accurate as the theoretical calculations say they should be.


Is there any way to do this?

Do you know the "Method of Manufactured Solutions"?

It is quite standard for code verification.

Regards

You test it by brute forcing it across several orders of magnitude and observing (i.e. back-calculating) the order of convergence. Some practical issues arise of course... it's really difficult to mesh one computational domain with mesh cells varying by several orders of magnitude and expecting to get "similar" results each time. That's why you prescribe as much as you can an analytical profile. That is, prescribe a velocity field and then use that to test only the discretization of the field. And that is precisely the basis of the method of manufactured solutions already mentioned.

I am not sure about the goal of your question ... Are you interested in checking the accuracy order of single part/subroutines of your code or the final global accuracy?

Being able to do both would be appreciated. I will look into Method of Manufactured Solutions, as that was recommended before.

I would first add that most of the times, a coding error is not going to show up as a different accuracy order but either crash your code or just give plainly wrong results (0th order or divergence). That is, it is typically by chance that, for certain tests, an error just gives a lower order, it is more probable that it in certain other tests the code more simply fails.

For the accuracy tests, the step 0, which is a special case of the method of manufactured solutions, is to just pick up know exact solutions to test the code against. The moms is step 1 and generally requires adding source terms to make an ideally arbitrary function a solution to your problem. However, if your code can handle non constant source terms this really is no different from step 0 (but typically requires you to work out a specific solution take derivatives, etc.)

The way both work is that you know that a certain analytical solution holds for a given set of boundary conditions (and possibly source terms for moms) then you just try to solve the given problem with given bc (and possibly source terms) with your code with different grid and/or time steps. As you know the analytical solution, you can then compute the error (which is a vector for each grid) and show that a given norm of this error reduces with respect to the grid and/or time step following the expected accuracy of the method you are using.

The more pieces of your code the analytical solution will test the more you gain confidence on the correctness of the code. You then typically put together as much of them as necessary by the code features.

This way of proceeding is ideal because you work at the whole code level. Still, sometimes you can't simply do that or it is not efficient. Then, at the cost of introducing specific unit tests, you simply apply the same methodology on smaller parts of your code. For example, if you have a runge kutta integration routine, you might want to test it out of your code with simpler ODEs for which you know the analytical solution.