%STREAMING STEP
for i=1:9 
    fIn(i,:,: ) = circshift(fOut(i,:,: ), [0,cx(i),cy(i)]); 
end

%STREAMING STEP
% implementation of the streaming without "periodicity"
fIn(2, 2:lx, : ) = fOut(2, 1:lx-1,: ); 
fIn(3, :, 2:ly) = fOut(3, :, 1:ly-1); 
fIn(4, 1:lx-1, :) = fOut(4,  2:lx , :); 
fIn(6, 2:lx, 2:ly) = fOut(6,  1:lx-1, 1:ly-1); 
fIn(5, :,  1:ly-1 ) = fOut(5,  :,  2:ly); 
fIn(7, 1:lx-1,  2:ly) = fOut(7, 2:lx , 1:ly-1 );
fIn(8, 1:lx-1, 1:ly-1) = fOut(8,  2:lx , 2:ly); 
fIn(9, 2:lx , 1:ly-1) = fOut(9,    1:lx-1, 2:ly);

relaxation paramters:
  s0=s3=s5=0.0;   s1=1.64;   s2=1.54;   s4=s6=1.7;   s7=s8=1.0/tau;

macroscopic value:
  rho=\sum f_i;
  j = \sum f_i*e_i + 0.5 * F;  // flux, updated with half force 

equilibrium values of the moments:
  m_eq0= rho; m_eq1= -2.0*rho+3.0*(j_x^2+j_y^2)/rho; ....

transfer from population space to moment space:
  m =  T f;  // T is the matrix

collsion in moment space:
  
 m*= s (m_eq-m);

j" = j'+0.5*F; // another half part of the force.

inverse transfer from moment space to population space:
  f *=T^{-1} m*

End of the collision step

bool steadyState2Dc ( MultiBlockLattice2D<T,DESCRIPTOR>& lattice, Box2D domain, T threshold, plint component, T& diff) {
    static MultiTensorField2D<T,2> uOld(domain.x1+1, domain.y1+1);
    MultiTensorField2D<T,2> uNew(domain.x1+1, domain.y1+1);
    computeVelocity(lattice, uNew, domain);
    subtractInPlace(uOld, uNew, domain);
    divideInPlace(uOld, uNew, domain);

    T diffmax = fabs(computeMax(*extractComponent(uOld, domain, component)));
    T diffmin = fabs(computeMin(*extractComponent(uOld, domain, component)));

    uOld = uNew;
    
    diff = (diffmax <= diffmin) ? diffmin : diffmax;

    return (diff < threshold);
}