Problemas al implementar el método Lattice Boltzmann

Quiero hacer una implementación 2D básica del método Lattice Boltzmann (en javascript) para simular el flujo de gas/fluido, pero me he encontrado con algunos problemas y parece que no puedo encontrar la causa.

El problema es que cuando ejecuto el paso de colisión (todavía no he implementado completamente la transmisión), la simulación explota y se ejecuta hasta el infinito y revisé mi código minuciosamente en busca de errores.

Comienzo con la inicialización de mis variables:

function Grid(width, height){
        this.width = width;
        this.height = height;
        this.f = [];
        this.e = [
            new vec2(0,0),
            new vec2(0,1),
            new vec2(1,0),
            new vec2(0,-1),
            new vec2(-1,0),
            new vec2(1,1),
            new vec2(1,-1),
            new vec2(-1,-1),
            new vec2(-1,1)
        ];
        this.w  = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36];
        for(var x=0; x<width; x++){
            this.f[x] = [];
            for(var y=0; y<height; y++){
                this.f[x][y] = [];
                for(var i=0; i<9; i++){
                    this.f[x][y][i] = 20+Math.random();
                }
            }
        }
    }

Aquí this.f es la función de distribución con 9 direcciones: 0 para ningún movimiento, 1-4 para movimiento recto y 5-9 para velocidades inclinadas. La función de distribución tiene valores que van desde 20-21. This.e contiene todas las direcciones como vectores (he hecho mi propia clase simple vec2) y this.w contiene los pesos para cada dirección. La primera pregunta que quiero hacer es si inicié la distribución con los valores correctos. No pude encontrar mucha información sobre los valores de inicialización.

En segundo lugar, mi (s) error (es) podría (n) estar en el paso de colisión:

    collide: function(){
        for(var x=0; x<width; x++){
                for(var y=0; y<height; y++){
                    var rho = 0;
                    for(var i=0; i<9; i++){
                        rho += this.f[x][y][i];
                    }

                    var u = new vec2(0,0);
                    for(var i=0; i<9; i++){
                        u  = u.add( this.e[i].scale( this.f[x][y][i] ) );
                    }
                    u = u.scale( 1/rho );
                    for(var i=0; i<9; i++){
                        var eu = u.dot(this.e[i]);
                        var uu = u.dot(u);
                        var fiEq = this.w[i]*rho*(1+3*eu+4.5*eu*eu-1.5*uu);
                        this.f[x][y][i] = this.f[x][y][i] - .005*(this.f[x][y][i]-fiEq);
                    }
                }
            }
        }

Primero calculo la densidad macroscópica rho, luego, con la suma de vectores calculo la velocidad macroscópica u. Los uso en la siguiente fórmula para calcular la distribución de equilibrio: y luego relajo la distribución hacia el equilibrio. como segunda pregunta me gustaría preguntar si implementé correctamente el paso de colisión.

Un ejemplo en vivo se puede ver aquí

Si necesita más información para ayudar, hágamelo saber y si pudiera ayudarme, ¡realmente lo agradecería!