
% Version date: 25. october 2013

%% Install MRST (Matlab Reservoir Simulator Toolbox)
% Follow the <http://www.sintef.no/Projectweb/MRST/Downloadable-Resources/Download/ instructions>
%

%% Tutorials
% Have a look to the tutorials at <http://www.sintef.no/Projectweb/MRST/Tutorials/ tutorials>

%
% Setup data
%

% Create a 2D cartesian grid.
nx = 40;
ny = 40;
G = cartGrid([nx ny 1]);
G = computeGeometry(G);

% We neglect gravity
gravity off;

% Number of faces
nf = G.faces.num;

% Number of faces
nc = G.cells.num;

perm_distribution = 'random'; % 2 cases are implemented: 'uniform' and 'random'

switch perm_distribution
  case 'uniform'
    rock.perm = repmat(100*milli*darcy, [G.cells.num, 1]); % "typical" permeability
  case 'random'
    % rock.perm = max(1, (10*ones(nc, 1) + 200* randn(nc, 1)))*milli*darcy;
    perm_data = load('perm.mat');
    rock.perm = perm_data.perm;
end
mu = 1*centi*poise; % "typical" viscosity
rock.perm = rock.perm/mu;  % "effective" permeability
rock.poro = repmat(0.3, [G.cells.num, 1]);

%
% Setup pressure solver and solve for pressure.
%

% Compute half transmissibilities
hT =  computeTrans(G, rock);
% and, from the half-transmissibilities, we compute the transmissibilities
fT = 1 ./ accumarray(G.cells.faces(:,1), 1 ./ hT, [G.faces.num, 1]);

% Identify internal faces
i  = all(G.faces.neighbors ~= 0, 2);

% We assemble the matrix T
% We have u_ij = -T_ij(p_j - p_i)
I  = [G.faces.neighbors(i,1); G.faces.neighbors(i,2)];
J  = [G.faces.neighbors(i,2); G.faces.neighbors(i,1)];
T  = sparse(double(I), double(J), [fT(i); fT(i)], nc, nc);

% Conservation of mass in each cell i gives nc equations:  sum_j u_ij = 0 (for nc unknown pressures)
d  = accumarray(reshape(G.faces.neighbors(i,:), [], 1), repmat(fT(i), [2,1]),  [nc, 1]);
A  = -T + sparse(1:nc, 1:nc, d, nc, nc);

% We set the pressure in the middle (injectors) and at the corners (producers)
% middle cell number
middle_cell = round(nx*(ny/2) + nx/2);
% corner cell number
corner_cells = [1; nx; nx*(ny - 1) + 1; nx*ny];

pcells = [middle_cell; corner_cells];
npcells = numel(pcells);
pvalues = [300; 200*ones(4, 1)]*barsa;

% We incorporate the given pressure in the matrix A.
A = A + sparse(pcells, pcells, ones(npcells, 1));
rhs = sparse(pcells, ones(npcells, 1), pvalues ,nc, 1);

% We clean up ...
clear I J V C D d;

% We compute the pressure.
p = A\rhs;



%
%  Visualization of results.
%


% We plot the permeability field
figure(1)
clf
plotCellData(G, rock.perm/(milli*darcy))
h = colorbar;
set(get(h,'ylabel'), 'string', 'Permeability in mD')
do_save_files = true;
if do_save_files
   print('-dpng', '-r300', sprintf('perm_%s.png', perm_distribution));
end


% We plot the pressure field
figure(2)
clf
plotCellData(G, p/barsa)
h = colorbar;
set(get(h,'ylabel'), 'string', 'Pressure in bar')
if do_save_files
   print('-dpng', '-r300', sprintf('press_%s.png', perm_distribution));
end


%
% Computation of streamlines and time of flight.
%


% We compute an approximation of outgoing flux in each cell.
ni   = G.faces.neighbors(i,:);
flux = -accumarray(find(i),  fT(i) .*(p(ni(:,2))-p(ni(:,1))), [nf, 1]);
x_outflux = accumarray(G.faces.neighbors(i,1), full(max(flux(i), 0).*sign(G.faces.normals(i,1))), [nc, 1]);
x_outflux = x_outflux + accumarray(G.faces.neighbors(i,2), full(min(flux(i), 0).*sign(G.faces.normals(i,1))), [nc, 1]);
y_outflux = accumarray(G.faces.neighbors(i,1), full(max(flux(i), 0).*sign(G.faces.normals(i,2))), [nc, 1]);
y_outflux = y_outflux + accumarray(G.faces.neighbors(i,2), full(min(flux(i), 0).*sign(G.faces.normals(i,2))), [nc, 1]);


% We plot the outgoing fluxes.
figure(3)
clf
x = G.cells.centroids(:,1);
y = G.cells.centroids(:,2);
quiver(x, y, x_outflux, y_outflux);

mrstModule add streamlines  % Add streamlines module from mrst.

% We set up for computation of time of flight.
p = full(p);
flux = full(flux);
state.flux = flux;

% Starting positions for the streamlines: We consider a kind of ring around the producer.
left_cells  = nx*(floor(ny/3):2*floor(ny/3)) + floor(nx/3);
right_cells = nx*(floor(ny/3):2*floor(ny/3)) + 2*floor(nx/3);
top_cells   = (floor(nx/3):(2*floor(nx/3))) + 2*nx*floor(ny/3);
bottom_cells   = (floor(nx/3):(2*floor(nx/3))) +  nx*floor(ny/3);

startpos = [left_cells(1:end-1), top_cells(1:end-1), right_cells(end:-1:2), bottom_cells(end:-1:2)]';

% We compute the time of flight for each streamline.

npos = numel(startpos);
plot_index = 1:5:npos; % We consider only a subset of all starting positions.
npi = numel(plot_index);


Tflight = zeros(npos, 1);
% forward streamline
[S_f, Tflight_f, C_f] = pollock(G, state, startpos);
% backward streamline
state.flux = -state.flux;
[S_b, Tflight_b, C_b] = pollock(G, state, startpos);
for i = 1:npos
   Tflight(i) = sum(Tflight_f{i}(Tflight_f{i} ~= Inf)) + sum(Tflight_b{i}(Tflight_b{i} ~= Inf));
end
Tflight = Tflight/day;
Tflight_viz = Tflight(plot_index);
if do_save_files
   save(sprintf('Tflight_%s.txt', perm_distribution), '-ascii', 'Tflight_viz');
end

% We compute the flux between the plotted streamlines.
outflux_between_streamlines = zeros(npi, 1);
for i = 1:npi
   ii = plot_index(i);
   if i < npi
      iip =  plot_index(i + 1) - 1;
   else
      iip = npos;
   end
   outflux_between_streamlines(i) = sum(abs(x_outflux(startpos(ii:iip)))) + sum(abs(y_outflux(startpos(ii:iip))));
end
outflux_between_streamlines = outflux_between_streamlines*day;
if do_save_files
   save(sprintf('obs_%s.txt', perm_distribution), '-ascii', 'outflux_between_streamlines');
end

%
%  Visalization of streamlines.
%

figure(3)
clf
hold on
% forward streamline
h1 = streamline({S_f{plot_index}});
% backward streamline
h2 = streamline({S_b{plot_index}});
h = [h1; h2];
set(h, 'LineWidth', 3);

for i = 1:npi
   ii = plot_index(i);
    text('position', [G.cells.centroids(startpos(ii),1), G.cells.centroids(startpos(ii),2), 3], ...
         'string',  sprintf('%d', i),  'fontsize', 16, ...
         'backgroundcolor', 'w')
    % text('position', [G.cells.centroids(startpos(ii),1), G.cells.centroids(startpos(ii),2), 10], ...
    %      'string',  sprintf('%d', i))
end

if do_save_files
   print('-dpng', '-r300',  sprintf('streamlines_%s.png', perm_distribution));
end




