¿Cómo se depura el código numérico? ¿Cuál podría ser la fuente de este error oscilatorio?

En silencio se puede obtener mucha información de la experiencia, me preguntaba si alguien ha visto algo similar a esto antes. El gráfico muestra la condición inicial (verde) para la ecuación de advección-difusión, luego la solución en la iteración 200 (azul) y luego nuevamente en la iteración 400 (rojo).

La solución de la ecuación de advección-difusión explota después de algunas iteraciones. El número de Péclet , y la condición CFL se cumple, , por lo que las ecuaciones deben ser estables. Anticipo que tengo un error en el código numérico.$\mu\approx0.07$$C\approx 0.0015$

Antecedentes. La discretización es la diferencia central para los términos de advección y difusión. Creo que este es el primer orden de advección y el segundo orden de difusión. He implementado esto utilizando un enfoque de volumen finito (por primera vez) en el que los valores de los coeficientes (velocidad y coeficiente de difusión) en las caras de la celda se encuentran por interpolación lineal de los promedios de la celda. Aplico la condición de límite de Robin en las superficies izquierda y derecha y establezco el flujo en los límites a cero.

¿Cómo depuras tu código numérico? ¿Alguien ha hecho una escena como esta antes, dónde sería un buen lugar para comenzar a buscar?

Actualizar

• Aquí están mis notas personales de estilo "libro de laboratorio" sobre la implementación de un método de volumen finito para la ecuación de advección-difusión, http://danieljfarrell.github.io/FVM/
• El código fuente de Python está disponible aquí, http://github.com/danieljfarrell/FVM.git

Actualizar

¡La solución no podría ser más simple! Acabo de hacer un error de señal en el término de difusión. Es extraño, estoy seguro de que no había publicado esto, ¡no habría encontrado el error! Si alguien quiere compartir consejos sobre cómo depurar su código numérico, todavía estoy interesado. No tengo un método, es un poco impredecible, sigo intentando cosas para obtener pistas, pero este proceso puede llevar semanas (a veces).

La prueba funciona ( tenga en cuenta que con el método de volumen finito todo lo que necesita hacer para calcular el área es una suma de ancho altura para todas las celdas, si utiliza un método de integración como numpy.trapz sus resultados incluyen el numérico error del método del trapecio). ¿Que está sucediendo aquí? Hay una velocidad constante y coeficientes de difusión pero con condiciones límite cerradas . Por lo tanto, en el límite vemos el equilibrio entre el campo de velocidad empujando hacia la derecha y el empuje de difusión hacia la izquierda.$\times$

He recopilado algo de mi experiencia en la depuración de códigos numéricos aquí: deal.II Preguntas frecuentes: depuración . No sé si eso te habría ayudado en este caso en particular, pero puede serlo en otros.