# If a numerical solution remains constant for different grid sizes, what does it mean?

I'm testing a finite volume scheme, Godunov type solver, using a problem with analytical solution. I'm not able to reproduce the solution, which includes source terms. I tried with different mesh sizes and I get the same result but it's not what I should get.

Because I suggested to one of my professors that maybe the code is not actually capable to reproduce this problem, as the paper says, so he asked me to try with different mesh sizes and he said that if it's a numerical thing, the numerical solution should be improved with the mesh refinement, but it didn't. How should I interpret those results?

This is the paper I'm working with: http://www.sciencedirect.com/science/article/pii/S0021999112003464

• There's a lot of loose ends in your question, here. But most signs point to simply a wrong implementation of the method. Have you tried your numerical method on a simpler problem? Could you give some more detail? Jun 14, 2017 at 4:24
• I'm solving the shallow water equations. I tried with two different test cases before, and I got the right results, but these didn't include some additional terms that this last case includes. Jun 14, 2017 at 16:13
• if the original implementation is fine, you should monitor the additional terms while solving the equation. maybe the differ by an order of magnitude? (by a bug or physics ) so you won't see the difference.
– Bort
Jul 6, 2017 at 16:10

Without more detail on what is "wrong" with your result I would say the follwing: The numerical implementation is not providing a "verification" result. In other words your implementation is not solving either the discretized equations correctly (e.g. a bug in your code) or the theoretical basis for your method is incorrect (e.g. the discretization is not consistant with the equations).

numerical solution should be improved with the mesh refinement, but it didn't. How should I interpret those results?

Unfortunately, all one can say is it means something's wrong.

Reducing $\Delta x$ should decrease the error (compared to the analytic solution). This is because techniques for numerically estimating derivatives get more accurate with smaller $\Delta x$. If your code was perfectly correct, then the only error would be this numerical error (i.e. error due to the grid and timestep not being infinitesimally fine). But because your result isn't improving with finer grids in this way, it shows there's "some other problem".

You can analytically verify that the given analytical solution really is a solution to the governing differential equations (i.e. substitute the solution into the governing equation, and work out the derivatives, and then you can confirm that you end up with those given source terms). A symbolic calculator is probably needed, like sympy (for python).

So, there's something wrong with your program - a bug. It could be anything, from the cfd code (e.g. how you discretize), to typos in your program.

Apart from standard debugging, code review etc, one way to narrow it down is to simplify the governing equation, by removing some terms. Then, substitute the same given analytic solution from before into the governing equation - which should yield different source terms. (Because the solution is the same as before, the initial and boundary conditions stay the same). Now, do as you did before, but with the new source terms: run the simulation, and compare the result with the analytic solution. Each time you refine the grid-size (i.e. reduce $\Delta x$), the error should decrease. By seeing which of the simplified governing equations behave coreectly or incorrectly, you can narrow down which terms are causing the problem, and thus which part of the code and/or program are buggy.

BTW This is the Method of Manufactured Solutions, though there's usually more care with how much the error decreases with decreasing $\Delta x$ - it should be $O( (\Delta x)^p )$, where $p$ is the order of the estimation technique (for example finite difference "forward Euler" is 1st order, $p=1$; central difference is 2nd order, $p=2$).

Because programs usually have bugs, it's unlikely that this behaviour is not caused by a bug. But, if the code and your program are 100% bug-free, then it could be a deeper issue with this scheme and/or this solver, for this particular problem. But it's almost certainly a bug. So look into that first.

• I'm solving the shallow water equations, I don't know if you are familiar with those, but they include two source terms: one due to the thrust and the other one due to friction. I tried two different test cases including only thrust and I was able to reproduce the same results than the paper, but when I run the last one (which includes friction and thrust), I don't get what I should get. I know for sure that the implementation of the thrust term is correct and the friction is so simple that I find hard to have an error on the implementation. Jun 14, 2017 at 16:32
• However, I tried to run the test case without the friction, and it doesn't converge to something, at some point it crashes...when I run the test case only with friction, I get a nice plot (which is not the solution, of course, but I'm able to get something). But as I said, after I tested the code with the thrust only (with the other test cases), I used the same code to add only the friction term, that's why I don't think there is a mistake on that part. Jun 14, 2017 at 16:34
• What kind of bugs could cause this behaviour in a numerical implementation? I'm using matlab... Jun 14, 2017 at 16:36
• @Jaydi_21 Could you add that info about the source terms to the question? Also, the paper whose results you want to reproduce - the details might help. Although the friction is simple, I assume you've looked over that part of the code already. Sometimes, bugs are revealed in already-tested code when it is exercised differently. That is, there was a bug in the code, but it got the right result for some tests anyway. If you've reviewed the code already, then start looking where you wouldn't expect bugs to be. e.g. It''s possible there's a typo in the paper for this test case. Jun 14, 2017 at 18:06
• I was able to get a better fit of my solution, but still it's not the exact one, even if I do mesh refinement. According to the paper, the solver gives an "accurate" solution, that means it's not able to reproduce the exact solution, then! Is it a rule that the solution should improve as I do mesh refinement? Because I thought that a numerical solver is just an approximation, that might be or not be equal to the exact solution, although we wish they were all able to reproduce the exact solution. Aug 31, 2017 at 19:26