Code covered by the BSD License  

Highlights from
immiscible LB

immiscible LB

by

 

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

Contact us