Problems implementing Lattice Boltzmann method I want to make a basic 2D implementation of the Lattice Boltzmann method (in javascript) to simulate gas/fluid flow, but I have run into some problems and I can't seem to find the cause of it.
The problem is that when I run the collision step (I haven't fully implemented streaming yet), the simulation blows up and runs off to infinity and I have checked my code very thoroughly for mistakes.
I start with initialising my 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();
                }
            }
        }
    }

Here this.f is the distribution function with 9 directions: 0 for no movement, 1-4 for straight movement and 5-9 for the slanted velocities. The ditribution function has values ranging from 20-21. This.e contains all the directions as vectors (I have made my own simple vec2 class) and this.w containts the weights for each direction. The first question that I want to ask is did I initialise the distribution with the correct values? I couldn't find much information on the initialisation values.
Secondly my mistake(s) could be in the collision step:
    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);
                    }
                }
            }
        }

First I calculate the macroscopic density rho, then, with vector addition I calculate the macroscopic velocity u. I use those in the following formula to calculate the equilibrium distribution: 
(source: nallatech.com)
and then I relax the distribution towards equilibrium.
as a second question I would like to ask if I implemented the collision step properly.
A live example can be seen here
If more information is needed to help let me know and if you could help me out I would really appreciate it!
 A: A couple of comments on your implementation:


*

*You initialize your distribution functions in an unconventional way; usually, we want to convert known macroscopic 'lattice' quantities like $\rho$ and $\vec{u}$ into $f_i$. An easy way to do this is using the equilibrium distribution you define above:
$$f_{i,eq}\left(\rho,\vec{u}\right)=w_i\rho\left[1+\frac{\vec{e}_i\cdot\vec{u}}{c_s^2}+\frac{1}{2}\frac{\left(\vec{e}_i\vec{e}_i-c_s^2\boldsymbol{I}\right):\vec{u}\vec{u}}{c_s^4}\right]$$
Let's make our life simple and choose $\rho=1$ and $\vec{u}=0$ then initialize simply with:
$$f_{i} = f_{i,eq}\left(1,0\right) = w_i$$
Note: 


*

*In the 'lattice' system the value of the macroscopic quantities can be chosen (more or less) arbitrarily to suit our need as long as the relevant dimensionless numbers (e.g. Reynolds number) are equivalent in the 'dimensional' system. 

*Make sure that velocity is chosen small enough to ensure the Mach number $\mathrm{Ma}=u/c_s\ll1$; usually this is limited to $u\approx0.1$. This is due to the small $\mathrm{Ma}$ expansion form of the equilibrium distribution which is to $O\left(\mathrm{Ma}^2\right)$ terms. My guess is your code diverges because the values chosen for $f_i$ in the initialization correspond to a velocity which is way too high.

*For initialy transient flows, this method of initializing is not the best way but for $\vec{u}\left(x,0\right)=0$ it is fine.


*Your implementation of the collision step looks correct but your relaxation 'frequency' $\omega=0.005$ seems very small too me (i.e. $\tau=200$). Likely, this works fine but start with $\omega=1$ first as this is the most stable value (it means in one 'lattice' timestep we relaxe towards the equilibrium distribution). A small note: var uu = u.dot(u); can be taken outside of the loop as it does not depend on i.

*I doubt the code will work without an implementation for the streaming step. If you initialize with the equilibrium distribution as above then after the first timestep, the distribution becomes $f_i\left(x,t+1\right)=f_i\left(x,t\right)$ and stays that way for all subsequent timesteps. This is actually a good check to see if the collision step was implemented correctly.

*I sense you have recently started working with LBM so let me give you some literature which deal with the numerical aspects:


*

*Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers, by M.C. Sukop and D.T. Thorne

*Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes, by A.A. Mohamad

*The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, by S. Succi



Good luck
