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

Red_and_Blue_Stream_D2Q9_2011.m
         
%%  Red=gdouble(Red); Blue=gdouble(Blue); % jacket test move to CPU
  
            tempR=Red; 
            tempB=Blue;

       i32_1=int32(1);
       
% Forward Propagation step &  Bounce Back (collision fluid with obstacles)

%sizBl=size(Blue);
sizBl=[Nr , Mc, 9];
%Blue(sizBl)=Blue(sizBl); % preallocation trick 
%Red(sizBl)=Red(sizBl); % preallocation trick 

for ic=1:1:N_c-1, % select velocity layer
    
    %ic=uint32(ic);
    
    ic2=ic_op(ic); % selects the layer of the velocity opposite to ic for BB
   
    NMic1= NxM*(ic-1); NMmNr= NMic1-Nr;
    % VECTORIAL
    % section 1: cares for wet locations that are NOT on the border 
  
                        % internal location i.e. away from borders 
                        % i=iawint; j=awint; % 
                        
                        i2 = iawint+C_y(ic); 
                        j2 = jawint+C_x(ic); % new i,j location 
                        i2 = yi2(i2+i32_1); %j2 = j2(j2+1); % new i,j locations corrected for PBC
                        
                        %i2j2a_ic = (j2-int32(1)).*Nr + i2 + NMic1; 
                        %i2j2a_ic = bsxfun(@times,(j2-1),Nr) + bsxfun(@plus,i2,NMic1); 
                        
                        i2j2a_ic0 = j2.*Nr+i2+NMmNr; 
                        idx=ijaint+NMic1; % ijaint defined in pre computed absolute indeces
                        
                       %%%%%% old %%%%
                            
                         %   Blue(sizBl)=Blue(sizBl); % preallocation trick
                         %  Blue(i2j2a_ic0)=tempB(idx);
                            
                          %  Red(sizBl)=Red(sizBl); % preallocation trick
                           % Red(i2j2a_ic0) =tempR(idx);
                            
                        %%%%%%%%%%%%%%%%%%%%%%%%%%    
                        
     % section 2: cares for wet location touching the border
                       
                        i2 = iobs+C_y(ic); j2 = jobs+C_x(ic);
                        i2=yi2(i2+1); %j2 = j2(j2+1); % Per BC
                        % i2j2a where to go i.e. new location 
                        
                        i2j2a = (j2-1).*Nr+ i2; % new can be wet or dry 
                        tru_Cha = (Channel2D(i2j2a)==0) ; % dry subset of new loc , logical 
                        
                        i2j2a_ic = i2j2a+NMic1 ; % new loc for state dep (same state)  
                        ija_ic2  = (ijobs + NxM*(ic2-1)); % old loc but new state 
                       
                        % from old location going into new dry , then
                        % change state only
      
                      %  Red([ija_ic2(tru_Cha);i2j2a_ic0])  = tempR([ijobs(tru_Cha)+NMic1;idx]); 
                      %  Blue([ija_ic2(tru_Cha);i2j2a_ic0]) = tempB([ijobs(tru_Cha)+NMic1;idx]); 
                       % Blue(ija_ic2(tru_Cha)) =
                       % tempB(ijobs(tru_Cha)+NMic1); % old 

                       
                        % from old location going into new wet , then same
                        % state but move
                        notruC=~tru_Cha;
                                             
                       % Blue(i2j2a_ic(notruC)) = tempB(ijobs(notruC)+NMic1); 
                       % Red(i2j2a_ic(notruC))  = tempR(ijobs(notruC)+NMic1); 
                       % Red(sizBl)=Red(sizBl); % preallocation trick
                       
                        Red([ija_ic2(tru_Cha);i2j2a_ic0;i2j2a_ic(notruC) ])  = tempR([ijobs(tru_Cha)+NMic1;idx;ijobs(notruC)+NMic1]); 
                        Blue([ija_ic2(tru_Cha);i2j2a_ic0; i2j2a_ic(notruC)]) = tempB([ijobs(tru_Cha)+NMic1;idx;ijobs(notruC)+NMic1]); 
                        
                  
end ; %  for ic direction

%%Red=double(Red); Blue=double(Blue); % jacket test bring back to CPU

Contact us