% CIV 2552 - Mét. Num. em Problemas de Fluxo e Transporte em Meios Porosos
% 2010.1
% Exemplo de Fluxo em Aquifero Confinado 2D
% Fluxo prescrito no bordo esquerdo e carga hidraulica prescrita no bordo
% direito.
% Carga hidraulica prescrita nos bordos inferior e superior.
% Metodo dos Volumes de Controle em Diferencas Finitas
% Tecnica de Direções Alternadas 

clc;
clear all;

% Parametros do arquifero
Lx = 800.0;             % comprimento do arquifero em x [m]
Ly = 400.0;             % comprimento do arquifero em y [m]
e = 20.0;               % espessura do arquifero [m]
k = 1.0e-2;             % coeficiente de permabilidade [m/s]
Ss = 2.0e-4;            % armazenamento específico [1/m]
alpha = k/Ss;           % parametro de difusividade [m2/s]
sd = 0.0;               % fonte externa distribuida [m/s]
ql = 0.0;               % BC: fluxo prescrito em x = 0 [m/s]
ur = 20.0;              % BC: carga hidraulica prescrita em x = Lx [m]
ub = 20.0;              % BC: carga hidraulica prescrita em y = 0 [m]
ut = 40.0;              % BC: carga hidraulica prescrita em y = Ly [m]
u0 = 20.0;              % condicao inicial [m]
umin = ub;              % carga hidraulica minima [m]
umax = ut;              % carga hidraulica maxima [m]

% Parametros de discretizacao
nc = 31;                % numero de colunas (dim. na direção x)
nr = ceil(nc*Ly/Lx);    % numero de linhas (dim. na direção y)
Dx = Lx/(nc-1);         % tamanho da celula na direcao x [m]
Dy = Ly/(nr-1);         % tamanho da celula na direcao y [m]
Dt = 2.0;               % passo de tempo [s]
rx = alpha*Dt/(Dx^2)    % fator rx
ry = alpha*Dt/(Dy^2)    % fator ry
fx = (Dx^2)/(Dy^2);     % fator fx
fy = 1.0/fx;            % fator fy
Nsteps = 1000;          % numero total de passos de tempo
NPrintSteps = 50;       % numero de passos para tracar graficos

% Dimensiona e inicializa vetores e matrizes
u = ones(nr,nc)*u0;     % carga hidraulica nas celulas
s = ones(nr,nc)*sd;     % fonte externa nas celulas
d = zeros(nc,nr);       % matriz dos valores RHS (right hand size) dos
                        % sistemas de equacoes tridiagonais
x = (0.0:1.0:nc-1)*Dx;  % vetor com as coordenadas x dos pontos
y = (0.0:1.0:nr-1)*Dy;  % vetor com as coordenadas y dos pontos

% Monta os vetores constantes dos sistemas tridiagonais para
% solucao na direcao x (linhas de celulas):
% Cada linha tem dimensao nc.
% "b" (bottom) valido para (i=1)
% "g" (generico) valido para (i=2:nr-1)
% "t" (top) valido para (i=nr)
axb = zeros(nc,1);
bxb = ones(nc,1);
cxb = zeros(nc,1);
axg = (-1.0)*ones(nc,1);
axg(1) = 0.0;   % nao usado
axg(nc) = 0.0;
bxg = (2.0 + 1.0/rx)*ones(nc,1);
bxg(nc) = 1.0;
cxg = (-1.0)*ones(nc,1);
cxg(1) = -2.0;
cxg(nc) = 0.0;  % nao usado
axt = zeros(nc,1);
bxt = ones(nc,1);
cxt = zeros(nc,1);

% Monta os vetores constantes dos sistemas tridiagonais para
% solucao na direcao y (colunas de celulas):
% Cada coluna tem dimensão nr.
% "l" (left) valido para (j=1)
% "g" (generico) valido para (j=2:nc-1)
% "r" (right) valido para (j=nc)
ayl = (-1.0)*ones(nr,1);
ayl(1) = 0.0;   % nao usado
ayl(nr) = 0.0;
byl = (2.0 + 1.0/ry)*ones(nr,1);
byl(1) = 1.0;
byl(nr) = 1.0;
cyl = (-1.0)*ones(nr,1);
cyl(1) = 0.0;
cyl(nr) = 0.0;  % nao usado
ayg = (-1.0)*ones(nr,1);
ayg(1) = 0.0;   % nao usado
ayg(nr) = 0.0;
byg = (2.0 + 1.0/ry)*ones(nr,1);
byg(1) = 1.0;
byg(nr) = 1.0;
cyg = (-1.0)*ones(nr,1);
cyg(1) = 0.0;
cyg(nr) = 0.0;  % nao usado
ayr = zeros(nr,1);
byr = ones(nr,1);
cyr = zeros(nr,1);

% Obtem tamanho da tela
screen_sizes = get(0,'ScreenSize');

% Plota grafico da condicao inicial na linha media
fig1 = figure;
fig1_pos = get(fig1,'Position');
fig1_pos(1) = 0;   % posiciona figura no bordo esquerdo da tela
set( fig1,'Position',fig1_pos);
plot(x,u(ceil(nr/2),:));
xlabel('x');
ylabel('h');
title('Fluxo em Aquifero Confinado 2D: carga hidraulica na linha media');
hold on

% Desenha contour e linhas de fluxo da condição inicial
fig2 = figure;
fig2_pos = get(fig2,'Position');
fig2_pos(1) = screen_sizes(3) - fig2_pos(3); % posiciona no bordo direito
set(fig2,'Position',fig2_pos);

% Prepara as variaveis para o contour
[X,Y] = meshgrid(x,y); % cria uma malha retangular
inter = [umin:(umax-umin)/20.0:umax];  % intervalos do contour

[c ch] = contourf(X,Y,u,inter);  % cria contour de faixas coloridas
colorbar;
set(ch,'edgecolor','y');
[vx,vy] = gradient(u,-Dx/k,-Dy/k);
linhas_fluxo = streamslice(X,Y,vx,vy);
set(linhas_fluxo,'color','k');
axis tight;

% Inicia o cálculo no tempo
for itime = 1:Nsteps

 % Solucao das linhas (direcao x)
 % Atualiza a matriz RHS referente a solucao das linhas

 % Primeira linha (i=1): carga hidraulica prescrita no bordo inferior
 d(1,:) = ub;

 % Da segunda a penultima linha (i=2:nr-1): linha generica
 aux1 = circshift(u,[+1 0]); % aux1 = u(i-1,:)
 aux2 = circshift(u,[-1 0]); % aux2 = u(i+1,:)
 % (j=1): fluxo prescrito no bordo esquerdo
 d(2:nr-1,1) = fx*aux1(2:nr-1,1) ...
             + (1.0/rx - 2.0*fx)*u(2:nr-1,1) ...
             + fx*aux2(2:nr-1,1) ...
             + s(2:nr-1,1)*(Dx^2)/(k*e) - ql*2.0*Dx/k;
 % (j=2:nc-1): coluna generica
 d(2:nr-1,2:nc-1) = fx*aux1(2:nr-1,2:nc-1) ...
                  + (1.0/rx - 2.0*fx)*u(2:nr-1,2:nc-1) ...
                  + fx*aux2(2:nr-1,2:nc-1) ...
                  + s(2:nr-1,2:nc-1)*(Dx^2)/(k*e);
 % (j=nc): carga hidraulica prescrita no bordo direito
 d(2:nr-1,nc) = ur;

 % Ultima linha (i=nr): carga hidraulica prescrita no bordo superior
 d(nr,:) = ut;

 % Resolve primeira linha
 u(1,:) = TriDiagonal( nc, axb, bxb, cxb, d(1,:) );
 for i = 2:nr-1
  % Resolve linha generica
  u(i,:) = TriDiagonal( nc, axg, bxg, cxg, d(i,:) );
 end 
 % Resolve ultima linha
 u(nr,:) = TriDiagonal( nc, axt, bxt, cxt, d(nr,:) );
 % Fim da solucao das linhas

 % Solucao das colunas (direcao y)
 % Atualiza a matriz RHS referente a solucao das colunas

 % Primeira coluna (j=1): fluxo prescrito no bordo esquerdo
 % (i=1): carga hidraulica prescrita no bordo inferior
 d(1,1) = ub;
 % (i=2:nr-1): linha generica
 d(2:nr-1,1) = (1.0/ry - 2.0*fy)*u(2:nr-1,1) ...
             + 2.0*fy*u(2:nr-1,2) ...
             + s(2:nr-1,1)*(Dy^2)/(k*e) - ql*fy*2.0*Dx/k;
 % (i=nr): carga hidraulica prescrita no bordo superior
 d(nr,1) = ut;

 % Da segunda a penultima coluna (j=2:nc-1): coluna generica
 aux1 = circshift(u,[0 +1]); % aux1 = u(:,j-1);
 aux2 = circshift(u,[0 -1]); % aux2 = u(:,j+1);
 % (i=1): carga hidraulica prescrita no bordo inferior
 d(1,2:nc-1) = ub;
 % (i=2:nr-1): linha generica
 d(2:nr-1,2:nc-1) = fy*aux1(2:nr-1,2:nc-1) ...
                  + (1.0/ry - 2.0*fy)*u(2:nr-1,2:nc-1) ...
                  + fy*aux2(2:nr-1,2:nc-1) ...
                  + s(2:nr-1,2:nc-1)*(Dy^2)/(k*e);
 % (i=nr): carga hidraulica prescrita no bordo superior
 d(nr,2:nc-1) = ut;

 % Ultima coluna (j=nc): carga hidraulica no bordo direito
 d(:,nc) = ur;

 % Resolve primeira coluna
 u(:,1) = TriDiagonal( nr, ayl, byl, cyl, d(:,1) );
 for j = 2:nc-1
  % Resolve coluna generica
  u(:,j) = TriDiagonal( nr, ayg, byg, cyg, d(:,j) );
 end
 % Resolve ultima coluna
 u(:,nc) = TriDiagonal( nr, ayr, byr, cyr, d(:,nc) );
 % Fim da solucao das colunas

 % Plota solucao na linha media somente nos passos de impressao
 if( mod(itime,Nsteps/NPrintSteps) == 0.0 )
  figure(fig1);
  tempo = itime*Dt
  plot(x,u(ceil(nr/2),:));
 end

 % Desenha contour e linhas de fluxo da solucao corrente
 if( mod(itime,Nsteps/NPrintSteps) == 0.0 )
  figure(fig2);
  [c ch] = contourf(X,Y,u,inter);
  colorbar;
  set(ch,'edgecolor','y');
  [vx,vy] = gradient(u,-Dx/k,-Dy/k);
  linhas_fluxo = streamslice(X,Y,vx,vy);
  set(linhas_fluxo,'color','k');
  axis tight;
 end
end
