clear all
close all
clc

set(0,'DefaultLineLineWidth',2);
set(0,'DefaultAxesFontSize',16);
set(0,'DefaultLineMarkerSize',10);

%%%%%%%%%%%%% Physical parameters
L = 3;
Phi_w = 0;   % Boundary condition
Gamma = 0.01;
rho = 1;
u = 0.1;
t_End = 10;
 
%%%%%%%%%%%%% Numerical parameters
nx = 501;

% Specify number of time steps OR size of time step
%nt = 900;    dt = t_End/nt,
 dt = 0.1,    nt = t_End/dt,
 
theta = 0.5;   % for ex. 0 (explicit), 0.5 (Crank-Nicolson) or 1 (implicit)

%%%%%%%%%%%%%% Grid generation
dx = L/(nx-1);
Dx = dx;

x0 = 0:dx:L;
time = 0:dt:t_End;

check_dt_max_diffusion  = rho*dx^2/Gamma,
check_dt_max_convection = dx/u,

%%%%%%%%%%%%%% Initial condition: Gaussian
Phi_init = 40 * exp(-(x0-0.5).^2/0.2^2); 

%%%%%%%%%%%%%% Initialize the matrix
A = zeros(nx,nx);
b = zeros(nx,1);
for i=2:nx-1
    A(i,i-1) =              -theta*Gamma/dx - theta*rho*u; % W
    A(i,i+1) =              -theta*Gamma/dx;               % E
    A(i,i)   = rho*dx/dt + 2*theta*Gamma/dx + theta*rho*u; % P
end

% Boundary conditions
A(1,1)   = 1;
b(1)     = Phi_w;
A(nx,nx) = 1;
b(nx)    = Phi_w;

%%%%%%%%%%%%%% Time marching
nt = length(time);
Phi = zeros(nx,nt);
Phi(:,1) = Phi_init;

col = hot(12);
j = 1;
figure('color','w'), plot(x0, Phi(:,1), 'color',col(j,:)),
hold on, grid on, box on,
xlabel('x'), ylabel('\phi'), ylim([0 40]), drawnow, %pause
for k=2:nt    
    for i=2:nx-1        
        b(i) = Phi(i+1,k-1) *               (1-theta)*Gamma/dx...
             + Phi(i-1,k-1) *              ((1-theta)*Gamma/dx +(1-theta)*rho*u)...
             + Phi(i,  k-1) * (rho*dx/dt -2*(1-theta)*Gamma/dx -(1-theta)*rho*u);
    end
        
    % Solve
    Phi(:,k) = A\b;  
    
    if (mod(k, round(nt/5))==0)
        %hold off, 
        j = j+1;
        plot(x0, Phi(:,k), 'color',col(j,:)), hold on, grid on, box on,
        xlabel('x'), ylabel('\phi'), ylim([0 40]), drawnow
    end
end
set(gca,'xtick',0:3)
if (Gamma==0)
    exportgraphics(gca,'convection-snapshots.pdf')
else
    exportgraphics(gca,'convection_diffusion-snapshots.pdf')
end

%%%%%%%%%%%%%% Plot solution T(x,t)
figure('color','w','position',[874   378   560   420]), hold on
xlabel('t'), ylabel('x'),
surf(time,x0, Phi, 'edgecolor','none')
c = colorbar;   c.Label.String = '\phi';
set(gca,'ytick',0:3)
drawnow
if (Gamma==0)
    exportgraphics(gca,'convection-space_time_evolution.pdf')
else
    exportgraphics(gca,'convection_diffusion-space_time_evolution.pdf')
end

%%%%%%%%%%%%%% Plot solution T(t) at x=xp
xp = 0.75;
[~,ind_x] = min(abs(x0-xp)) % find index of node closest to x=xp
figure('color','w'), hold on, grid on, box on
plot(time, Phi(ind_x,:), 'b.')

ylim([0 40])
xlabel('t'), ylabel(['\phi(x=' num2str(x0(ind_x)) ')'])
if (Gamma==0)
    exportgraphics(gca,'convection-time_evolution.pdf')
else
    exportgraphics(gca,'convection_diffusion-time_evolution.pdf')
end

title(['L=1, Phi_w=' num2str(Phi_w)...
    ', \Gamma=' num2str(Gamma) ', \rho=' num2str(rho) ', n=' num2str(nx)...
    ', \Deltat=' num2str(dt) ', \theta=' num2str(theta)])