Code covered by the BSD License

# immiscible LB

by

### G. Sken (view profile)

23 Jul 2009 (Updated )

Implements Immiscible Lattice Boltzmann (ILB, D2Q9) method for two phase flows

Fast_Red_and_Blue_Collision_step.m
``` %Collision (between fluid elements) omega=1 , relaxation frequency

if(0) % with memory to waste without comp time advantages
rhoij=rho(ija); rhoijR=rhoR(ija);  rhoijB=rhoB(ija);

feq=bsxfun(@times,feq(:,(1:9)),(rhoij.^-1));

Red(all_indeces)= omegaR.*( feq(:).*repmat(rhoijR,9,1)) + (1-omegaR).*Red(all_indeces);

Blue(all_indeces)=omegaB.*( feq(:).*repmat(rhoijB,9,1)) + (1-omegaB).*Blue(all_indeces);

end

if(1) % with memory to waste without comp time advantages

rhoij=rho(ija); rhoijR=rhoR(ija);  rhoijB=rhoB(ija);
feq=bsxfun(@times,feq,(rhoij.^-1));
u_omR=1.-omegaR;
u_omB=1.-omegaB;

for ic=1:N_c; %
indeces=ija+NxM.*(ic-1);
%feq_ic=feq(:,ic);
Red([Nr , Mc, 9]) = Red([Nr , Mc, 9]); % pre-allocation trick
Red(indeces) = omegaR.*( feq(:,ic).*rhoijR) + u_omR.*Red(indeces);
Blue(indeces)= omegaB.*( feq(:,ic).*rhoijB) + u_omB.*Blue(indeces);
end

end

return

if(0) % works selectively in Red & Blue ... slower then the above

iR= rhoR(ija)>0; iB= rhoB(ija)>0 ; ijaR= ija(iR); ijaB=ija(iB);
rhoij=rho(ija); rhoijR=rhoR(ijaR); rhoijB=rhoB(ijaB);

feq=bsxfun(@times,feq,(rhoij.^-1));
u_omR=1.-omegaR; u_omB=1.-omegaB;

ii=(1:lena)'; ijR=ii(iR); ijB=ii(iB);

indecesR=ijaR;
indecesB=ijaB;

for ic=1:N_c; %

Red([Nr , Mc, 9])=Red([Nr , Mc, 9]); % pre-allocation trick
Red(indecesR) = omegaR.*( feq(ijR+lena.*(ic-1)).*rhoijR) + u_omR.*Red(indecesR);
Blue(indecesB)= omegaB.*( feq(ijB+lena.*(ic-1)).*rhoijB) + u_omB.*Blue(indecesB);
indecesR=indecesR+NxM;
indecesB=indecesB+NxM;
end

%%%Red=double(Red); Blue=double(Blue); % jacket test bring back to CPU
%%%feq=double(feq); rhoR=double(rhoR); rhoB=double(rhoB);%
end

return

if(0)
rhoij=rho(ija);rhoijR=rhoR(ija);rhoijB=rhoB(ija);

for ic=1:N_c; %-1
indeces=ija+NxM*(ic-1); indeces=uint32(indeces);
feq(:,ic) = feq(:,ic)./rhoij;
Red(indeces)= omegaR.*( feq(:,ic).*rhoijR) + ...
(1-omegaR).*Red(indeces);
Blue(indeces)=omegaB.*( feq(:,ic).*rhoijB) + ...
(1-omegaB).*Blue(indeces);
end

end

if(0)
rhoij=rho(ija);rhoijR=rhoR(ija);rhoijB=rhoB(ija);

for ic=1:N_c; %-1
indeces=ija+NxM*(ic-1); indeces=uint32(indeces);
f_indeces = f(indeces)./rhoij;
Red(indeces)= omegaR.*( f_indeces.*rhoijR) + ...
(1-omegaR).*Red(indeces);
Blue(indeces)=omegaB.*( f_indeces.*rhoijB) + ...
(1-omegaB).*Blue(indeces);
end

end

if(0) % slow !!!!

if Cur_Iter==2
%A=(repmat(int32(1:9) ,int32([lena,1])));
ind_ics= repmat(ija,int32([1,9]) )+NxM.*(repmat(int32(0:8) ,int32([lena,1])));
ind_ics=ind_ics(:); ind_ics=uint64(ind_ics);
ija_s= repmat(ija,[9,1]);
end

rhoij=rho(ija_s); rhoijR=rhoR(ija_s);  rhoijB=rhoB(ija_s);
f_indeces = f(ind_ics).*(rhoij.^-1);
Red(ind_ics)= omegaR.*( f_indeces.*rhoijR) + (1-omegaR).*Red(ind_ics);
Blue(ind_ics)=omegaB.*( f_indeces.*rhoijB) + (1-omegaB).*Blue(ind_ics);

end

%  omega==1
%  Red(ija+NxM*(ic-1))= f(ija+NxM*(ic-1)).*rhoR(ija)./rho(ija);
%  Blue(ija+NxM*(ic-1))=f(ija+NxM*(ic-1)).*rhoB(ija)./rho(ija);

%   if(0)  % slow
%   rhoij=rho;
%   f= f./repmat(rhoij,[1,1,9]);
%   Red= omegaR.*( f.*repmat(rhoR,[1,1,9])) + (1-omegaR).*Red;
%   Blue=omegaB.*( f.*repmat(rhoB,[1,1,9])) + (1-omegaB).*Blue;
%  end```