function grains_open_right_boundary (grain_arrangement_file, Nx, Ny, nu, rho_left, average_density, t_max)

% grains_open_right_boundary.m
%
% Created by Jeff Cunningham while on sabbatical from USF and visiting 
% CSIRO in Floreat, Western Australia, 2013-2014
%
% Uses lattice Boltzmann modeling to simulate 2D fluid flow through a
% porous medium.  The user specifies the geometry of the porous medium by
% providing the file grain_arrangement_file.
%
% Flow is left to right.
% This code differes from its predecessors, grains.m and 
%   grains_constant_rho_left_right.m, in the way it handles boundary
%   conditions.  This code uses a specified density (pressure) on the left
%   boundary, and uses an open right boundary.  The paper by Lou et al.
%   (2013, Physical Review E) discusses different implementations of the
%   open right bounda
%
% Some parts of this lattice Boltzmann code are modeled after codes posted
% on-line by Jonas Latt and co-workers.  I also have received assistance
% from Shadab Anwar and Mike Sukop as I have developed this code.  
%
% grain_arrangement_file is the name of a file containing data that specify
%   the geometry of the porous domain.  A 1 indicates that the node is
%   occupied by a grain; a 0 indicates that the node is fluid.
% Nx and Ny are the size of the domain in lattice units.  Nx and Ny should
%   agree with the data in grain_arrangement_file.
% nu is the fluid viscosity in LBM units (lu^2/ts).  
%   Recommend nu = 1/6.
% rho_left is the density on the left boundary.  Recommend a value of just
%   a little over 1.0, perhaps 1.005 or 1.01.  If rho_left is too high then
%   there could be issues with compressibility, i.e., mass balance.
% average_density is used to initialize the distribution functions.  It is
%   the density everywhere in the domain at the beginning of the simulations.
%   Recommend a value of 1.0.  Notice that average_density should probably
%   be a little lower than rho_left.
% t_max is the number of time steps for the simulation to run.
%   t_max must be high enough that the system reaches steady state.
%   As a rough guideline, use t_max = Nx*Ny/10 .




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Geometry of the problem    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid = fopen(grain_arrangement_file, 'r');
pore_wall = fscanf(fid, '%i', [Nx, Ny]);
fclose(fid); 
is_fluid = 1 - pore_wall;       % if it's not a pore wall, then it is fluid
[y, x] = meshgrid(1:Ny, 1:Nx);      

figure(1);
%imagesc(pore_wall');
surf(pore_wall');   shading flat;
axis('equal');
axis([1 Nx 1 Ny]);
xlabel('x coordinate (lattice units)', 'FontSize', 16);
ylabel('y coordinate (lattice units)', 'FontSize', 16);
title('Map of the porous medium', 'FontSize', 16);
drawnow;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Set some LBM variables    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
omega  = 1. / (3*nu + 1./2.);   % relaxation parameter
tau = 1/omega;

green_light = true;                % parameter to tell if we should proceed


    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    D2Q9 lattice constants    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Use this numbering scheme for the directions:
% 6 2 5
% 3 9 1
% 7 4 8
%
w  =  [1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36, 4/9];
opp = [  3,   4,   1,   2,    7,    8,    5,    6,   9];
%e_x = [ 1,   0,  -1,   0,    1,   -1,   -1,    1,   0];
%e_y = [ 0,   1,   0,  -1,    1,    1,   -1,   -1,   0];




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Pre-allocate the f matrices                       %%%%
%%%%    This makes Matlab happy                           %%%%
%%%%    Some of the matrices don't need pre-allocation    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = zeros(9, Nx, Ny);
f_temp = zeros(9, Nx, Ny);
%f_eq = zeros(9, Nx, Ny);   
%u_x = zeros(Nx, Ny);
%u_y = zeros(Nx, Ny);
%rho = zeros(Nx, Ny);

 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Initialize the distribution function f    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:9
    f(k, 1:Nx, 1:Ny) = w(k) * average_density;
end
initial_fluid_mass = sum(sum(   is_fluid .* reshape(sum(f),Nx,Ny)   ));    % Track the fluid mass over time.
disp(initial_fluid_mass);                                                   % Verifies code is running properly.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Begin main loop    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t = 0;

while (green_light) && (t < t_max)

    t = t+1;
       
   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%    Calculate macroscopic variables    %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [rho, u_x, u_y] = calc_macroscopic_variables_5 ( Nx, Ny, f );

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%    Every 50 time steps, check the mass balance    %%%%
    %%%%    Also display the velocity field                %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if (mod(t, 50) == 0)
        green_light = check_mass_balance(t, initial_fluid_mass, rho, is_fluid);
        u = reshape(sqrt(u_x.^2 + u_y.^2), Nx, Ny);  
        figure(2);  
        %imagesc(u');    
        %surf(u');   colorbar;    
        pcolor(u');
        axis('equal');    axis([1 Nx 1 Ny]);    shading interp;
        drawnow;        
    end  
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    %%%%    Calculate equilibrium distribution function     %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Use accel = 0, i.e., no external acceleration applied
    f_eq = calc_f_eq ( Nx, Ny, is_fluid, u_x, u_y, rho, 0, tau);
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%    Collision or relaxation    %%%%
    %%%%    Use standard BGK method    %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i = 1:Nx
        for j = 1:Ny
            if is_fluid(i,j)    % fluid nodes relax to equilibrium
                for k = 1:9
                    f_temp(k,i,j) = f(k,i,j) - (f(k,i,j) - f_eq(k,i,j))/tau;
                end
            else                % solid nodes use bounce-back
                for k = 1:9
                    f_temp(k,i,j) = f(opp(k),i,j); 
                end
            end
        end
    end
    

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%    Streaming step                                             %%%%
    %%%%    Top and bottom boundary conditions are periodic.           %%%%
    %%%%    Left boundary is specified density (pressure).             %%%%
    %%%%    Right boundary is open (Lou et al., 2013, Phys Rev E).     %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    f = stream_2(Nx, Ny, f, f_temp);
    f = periodic_top_bottom_2(Nx, Ny, f, f_temp);
    f = constant_rho_left_open_right_6 (Nx, Ny, f, u_x, rho_left);
    % I think it is important to put periodic_top_bottom_2 before
    % constant_rho_left_open_right_6, because of the way that the
    % constant-rho left boundary is handled.

end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Compute streamlines    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u_x = u_x';                         % Need to transpose these to get the 
u_y = u_y';                         % figure to turn out properly

% I have hard-wired the starting locations for the streamlines.  That might
% not be a great way to do this, but I have found that it seems to work
% pretty well for most possible grain configurations.
y_start = 3:6:Ny;
vec_length = length(y_start);
x_start = ones(vec_length, 1);
x_start = [x_start', 6:6:Nx, 6:6:Nx];
new_length = length(6:6:Nx);
y_start = [y_start, ones(1,new_length), (Ny-1)*ones(1,new_length)];
XY = stream2(x', y', u_x, u_y, x_start, y_start);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Graph streamlines    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3);
orient(3, 'portrait');
set(3, 'PaperPosition', [1 2.5 6.5 6]);
h = streamline(XY);
set(h, 'Color', 'green');
set(h, 'LineWidth', 1.);
axis('equal');
axis([1 Nx 1 Ny]);
xlabel('x coordinate (lattice units)', 'FontSize', 16);
ylabel('y coordinate (lattice units)', 'FontSize', 16);
title('Streamlines of flow through a porous medium', 'FontSize', 16);
hold on;
contour(x, y, pore_wall);
hold off;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Graph velocity vectors    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(4);
orient(4, 'portrait');
contour(x, y, pore_wall);
hold on;
h = quiver(x(4:4:Nx,4:4:Ny), y(4:4:Nx,4:4:Ny), (u_x(4:4:Nx,4:4:Ny))', (u_y(4:4:Nx,4:4:Ny))' );
set(h, 'Color', 'green');
xlabel('x coordinate (lattice units)', 'FontSize', 16);
ylabel('y coordinate (lattice units)', 'FontSize', 16);
title('Visualization of velocity vectors', 'FontSize', 16);
axis('equal');
axis([1 Nx 1 Ny]);
hold off;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    End the program    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end