%By liu-2017.0403. 谢谢这位没写名字的大佬。

%又上网扒代码了。基础太差，不看看大佬们的杰作是要把小的愁死啊~

%一样感人的效果。然而还不清楚这个是撒。。等我啃一啃再更新哈

%还是老样子，有啥问题Feel free to tell us~毕竟群众力量大嘛~QQ群：293267908。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
clear
% GENERAL FLOW CONSTANTS
lx         = 250;
ly         = 250;
obst_x = lx/5+1;   % position of the cylinder; (exact
obst_y = ly/2+1;   % y-symmetry is avoided)
obst_r = ly/10+1;  % radius of the cylinder
uMax  = 0.06;      % maximum velocity of Poiseuille inflow
Re     = 100;      % Reynolds number
nu    = uMax * 2.*obst_r / Re;   % kinematic viscosity
omega  = 1. / (3*nu+1./2.);      % relaxation parameter
maxT   = 4000;   % total number of iterations
tPlot  = 5;        % cycles

% D2Q9 LATTICE CONSTANTS
t  = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
cx = [  0,   1,  0, -1,  0,    1,  -1,  -1,   1];
cy = [  0,   0,  1,  0, -1,    1,   1,  -1,  -1];
opp = [ 1,   4,  5,  2,  3,    8,   9,   6,   7];
col = [2:(ly-1)];

[y,x] = meshgrid(1:ly,1:lx);
% obst = (x-obst_x).^2 + (y-obst_y).^2 <= obst_r.^2;
obst = (abs(y-4.*(x-110)-50)<=1.8).*(x>=110).*(x<=120)+...
       (abs(y+4.*(x-140)-50)<=1.8).*(x>=130).*(x<=140);
% obst = zeros(lx,ly);
obst([1,lx],:) = 1;
obst(:,ly) =1 ;
obst([1:110],1)=1;
obst([140:250],1)=1;
obst(110,[1:50])=1;
obst(140,[1:50])=1;
bbRegion = find(obst);

% INITIAL CONDITION: (rho=0, u=0) ==> fIn(i) = t(i)
fIn = reshape( t' * ones(1,lx*ly), 9, lx, ly);

% MAIN LOOP (TIME CYCLES)
for cycle = 1:maxT

    % MACROSCOPIC VARIABLES
    rho = sum(fIn);
    ux  = reshape ( ...
          (cx * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho;
    uy  = reshape ( ...
          (cy * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho;
   
    % MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS
      % Inlet: Poiseuille profile
    L = ly-2; y = col-1.5;
    ux(:,[111:139],1) = 0;
    uy(:,[111:139],1) = uMax;
    rho(:,[111:139],1) = 1 ./ (1-uy(:,[111:139],1)) .* ( ...
        sum(fIn([1,4,2],[111:139],1)) + ...
        2*sum(fIn([5,8,9],[111:139],1)));
%       % Outlet: Zero gradient on rho/ux
%     rho(:,lx,col) = rho(:,lx-1,col);
%     uy(:,lx,col)  = 0;
%     ux(:,lx,col)  = ux(:,lx-1,col);

    % COLLISION STEP
    for i=1:9
       cu = 3*(cx(i)*ux+cy(i)*uy);
       fEq(i,:,:)  = rho .* t(i) .* ...
         ( 1 + cu + 1/2*(cu.*cu) ...
           - 3/2*(ux.^2+uy.^2) );
       fOut(i,:,:) = fIn(i,:,:) - ...
         omega .* (fIn(i,:,:)-fEq(i,:,:));
    end

    % MICROSCOPIC BOUNDARY CONDITIONS
    for i=1:9
         % Left boundary
         fOut(i,1,col) = fEq(i,1,col) + ...
           18*t(i)*cx(i)*cy(i)* ( fIn(7,1,col) - ...
           fIn(6,1,col)-fEq(7,1,col)+fEq(6,1,col) );
%          % Right boundary
%          fOut(i,lx,col) = fEq(i,lx,col) + ...
%            18*t(i)*cx(i)*cy(i)* ( fIn(6,lx,col) - ...
%            fIn(9,lx,col)-fEq(6,lx,col)+fEq(9,lx,col) );
         % Bounce back region
         fOut(i,bbRegion) = fIn(opp(i),bbRegion);
    end

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

    % VISUALIZATION
    if (mod(cycle,tPlot)==0)
        u = reshape(sqrt(ux.^2+uy.^2),lx,ly);
        u(bbRegion) = nan;
        imagesc(u');
        axis equal off; drawnow
    end
end