%-----------------------------------------------------------------
% Script 'Discretization Methods', Example 3.5
% (c) 08/2004 by A. Schmidt, IAM, Uni Stuttgart
%-----------------------------------------------------------------
clear all

% Variables
l = 0.5;                            % length of the duct
T0 = 100.0;                         % bound. cond.: rel. temperature T(0) 
                                    % [T(l)=0]
v_a = 2.0e+4;                       % flow velocity v / diffusivity alpha 
m2 = 6.0e+4;                        % coefficient m^2
n = 51;                             % number of nodes (without ghost points)
n_a = 200;                          % number of sample points 
                                    % (analytical solution)

% Parameters
dx = l/(n-1.0);                     % nodal distance
Pe = v_a*dx;                        % Peclet number

%-----------------------------------------------------------------
% Calculation of the System Matrix A: Az=b [Eqs. (3.78), (3.80) or (3.83)]
% Coefficients of the central difference method (cd)
a_cd = 1.0+Pe/2.0; 
b_cd = -(2.0+m2*dx^2);
c_cd = 1.0-Pe/2.0;

% Coefficients of the upwind difference scheme (ud)
a_ud = 1.0+Pe; 
b_ud = -(2.0+m2*dx^2+Pe);
c_ud = 1.0;

%----------------------------
% central difference method
% first row (i=2)
A(1,1) = b_cd;
A(1,2) = c_cd;

% middle block
for i=3:1:n-1
  A(i-1,i-2) = a_cd;
  A(i-1,i-1) = b_cd;
  A(i-1,i) = c_cd;
end

% last row (i=n)
A(n-1,n-2) = a_cd;
A(n-1,n-1) = b_cd;

% Calculation of the right-hand side b: Az=b
b(1) = -a_cd*T0;
b(2:n-1) = 0.0;
b = b';                             % make b a column vector

%----------------------------
% Calculation of the vector of unknowns z (temperature): Az=b
z_cd = A\b;

%----------------------------
% upwind difference method
% first row (i=2)
A(1,1) = b_ud;
A(1,2) = c_ud;

% middle block
for i=3:1:n-1
  A(i-1,i-2) = a_ud;
  A(i-1,i-1) = b_ud;
  A(i-1,i) = c_ud;
end

% last row (i=n)
A(n-1,n-2) = a_ud;
A(n-1,n-1) = b_ud;

% Calculation of the right-hand side b: Az=b
b(1) = -a_ud*T0;
b(2:n-1) = 0.0;                      % (b is already a column vector)

%----------------------------
% Calculation of the vector of unknowns z (temperature): Az=b
% (vector b stays the same)
z_ud = A\b;

%-----------------------------------------------------------------
% Analytical solution
gamma2 = v_a/2.0-sqrt(m2+(v_a/2.0)^2);
for i=1:1:n_a
  x_a(i) = (i-1)*l/(n_a-1);         % sample point location
  z_a(i) = T0*exp(gamma2*x_a(i));
end

%-----------------------------------------------------------------
% Output
% sample points
for i=1:1:n
  x(i) = (i-1)*dx;
end

for i=n:-1:2                        % first sample point has to be added
  z_cd(i) = z_cd(i-1);
  z_ud(i) = z_ud(i-1);
end
z_cd(1) = T0;
z_ud(1) = T0;

% plot
figure(1)
plot(x,z_cd,'r',x,z_ud,'b',x_a,z_a,'k')
xlabel('Location / m')
ylabel('Temperature / °C')
title('red: cent diff      blue: upwind      black: analytical')