# How do you debug numerical code, what could be source of this oscillatory error?

Quiet a lot of insight can be gained form experience, I was just wondering if anybody has seen something similar to this before. The plot shows the initial condition (green) for the advection-diffusion equation, then the solution at iteration 200 (blue) and then again at iteration 400 (red).

The solution of the advection-diffusion equation blows up after a few iterations. The Péclet number $\mu\approx0.07$, and the CFL condition is satisfied, $C\approx 0.0015$, so the equations should be stable. I anticipate I have a bug in the numerical code.

Background. The discretisation is central difference for both advection and diffusion terms. I believe this is first order of advection and second order for diffusion. I have implemented this using a finite-volume approach (for the first time) in which the coefficients (velocity and diffusion coefficient) values at the cell faces is found by linear interpolation from the cell averages. I apply Robin boundary condition on the left and right surfaces and set the flux at the boundaries to zero.

How do you debug your numerical code? Has anybody scene something like this before, where would be a good place to starting looking?

Update

Update

The solution couldn't be more simple! I just made a sign error on the diffusion term. It's strange, I'm sure it I had not of posted this I would not have found the error! If someone wants to share tips on how they debug their numerical code I am still interested. I don't have a method, it's a bit hit and miss, I keep on trying stuff to get clues, but this process can take weeks (sometimes).

Proof it works (N.B. that with the finite-volume method all you need to do to calculate the area is a summation of width $\times$ height for all the cells, if you use an integration method such as numpy.trapz your results include the numerical error of the trapezium method). What is happening here? There is a constant velocity and diffusion coefficients but with closed boundary conditions. Therefore at the boundary we see the equilibrium between the velocity field pushing to the right and the diffusion push to the left.

• What sort of discretization are you using? What order method? What are your boundary conditions? – Geoff Oxberry Jun 21 '13 at 6:30
• Thank you @GeoffOxberry, I have updated with more details. Although to really understand what I have done you can read my lab book notes on the github link above. – boyfarrell Jun 21 '13 at 7:44
• I admire your conscentious approach to this problem. Are these plots from the time-dependent equations? If yes, what do you observe for $\theta = 0,0.5,1$? – Jan Jun 21 '13 at 8:54
• I tried that, it gives broadly the same. But now I feel stupid, I was testing for a sign error... setting $a=0$ and $d<0$ makes it work (in the diffusion limit, yes negative diffusion makes it work!). So for future reference, the above could be consistent with a sign error. Not completely tested yet so I will report back soon! PS. thanks for the compliment :) – boyfarrell Jun 21 '13 at 9:01