%  Finite differences for the advection-diffusion equation
%
%      - D u_xx   + v u_x   =   1
%
%     on the domain x=[0,1],
%     with boundary conditions u(0)=0, u(1)=0.
%
%     Computational Science and Engineering
%     Chapter 1, Section 1.2, p 25, Question 19.
%

%clear
%clc
%close all

% D=1; v=1; % try experimenting with these parameters
D=0.05; v=0.1;
% D=0.05; v=1; % boundary layer

Nmin=0; Nmax=4; 
Ns=logspace(Nmin,Nmax,round(Nmax-Nmin+1)*2); Ns=ceil(Ns); % must be integers
nN=numel(Ns);
hs=zeros(1,nN);
errorsU=zeros(1,nN); % upwind differences errors
errorsu=zeros(1,nN); % centered differences errors
for i=1:nN
    N=Ns(i)+1;
    xL=0; xR=1; h=(xR-xL)/(N-1); hs(i)=h; % h is the spatial discretization size
    x=xL:h:xR; x=x';
    
    % True solution:
    r=v/D;
    C=1/(v*(exp(r)-1));
    uTrue=x/v+C*(1-exp(r*x));
    
    % diffusion matrix
    e = ones(N,1); K = spdiags([-e 2*e -e], -1:1, N, N); K=D*K/h^2;
    %K=full(K)
    
    % centered difference matrix
    Del0 = spdiags([-e e], [-1 1], N, N); Del0=v*Del0/(2*h);
    %Del0=full(Del0)
    
    % upwind difference matrix
    DelPlus = spdiags([-e e], 0:1, N, N); DelPlus=v*DelPlus/h;
    %DelPlus=full(DelPlus)
    
    A=Del0+K; % centered advection matrix, plus diffusion matrix
    A(1,:)=0; A(end,:)=0; A(1,1)=1; A(end,end)=1; % boundary conditions
    b=ones(size(x));  b(1)=0; b(end)=0; % boundary conditions
    u=A\b;    % solve with backslash (centered)
    
    A=DelPlus+K; % upwind advection matrix, plus diffusion matrix
    A(1,:)=0; A(end,:)=0; A(1,1)=1; A(end,end)=1; % boundary conditions
    U=A\b;    % solve with backslash (upwind)
    
    errorsU(i)=norm(U-uTrue,inf); % upwind differences errors
    errorsu(i)=norm(u-uTrue,inf); % centered differences errors
end


% plot results
FS='FontSize';
figure(1)
lh=0;
loglog(hs, errorsu, 's', hs, errorsU, 'o', hs, exp(1*log(hs)+lh), '--', hs, exp(2*log(hs)+lh) )
xlabel('Spatial discretization size, h', FS, 16)
ylabel('Error=||u_{true} - u_{approx}||_{\infty}', 'FontSize', 16)
title({['- D u_{xx}   + v u_x   =   1 '],  ['   D=', num2str(D), '   v=', num2str(v)]}, FS, 16)
legend('u   centered', 'U   upwind', 'ref slope 1', 'ref slope 2')


disp('Finished.')

