% -------------------------------------------------------------------------------------- %
% Methodes de discretisation en fluides                                                  %
% Laboratoire II - Code de départ                                                        %
% M.A. Habisreutinger                                                                    %
% Jeudi 28 mars 2024                                                                     %
% -------------------------------------------------------------------------------------- %

function Lab02()
  
    % 1.   Parameters

    %  1.1  Mesh
    npoints = 16 ;
    limits = [0, +1] ;

    %  1.2  Diffusion and advection coefficients
    nu = 1e-02 ; 
    c = 1e+00 ;

    
    %  2.   Computation

    %  2.1  Mesh nodes
    x = linspace(limits(1), limits(2), npoints)' ;

    %  2.2  Mesh size
    h = (limits(2)-limits(1))/(npoints-1) ;

    %  2.3  Local Peclet number
    Pel = c*h/(2*nu) ;

    %  2.4  Analytical solution
    u = (exp(c*(x-1)/nu)-exp(-c/nu))/(1-exp(-c/nu)) ;

    %  2.5  Mass matrix (colocation method)
    e = ones(npoints, 1) ;
    M = spdiags(e, 0, npoints, npoints) ;

    %  2.6  Differentiation matrices (central finite differences - 2nd order)
    gradient = (1/(2*h))*spdiags([-1*e, 0*e, +1*e], [-1, 0, +1], npoints, npoints) ;
    laplacian = (1/(h^2))*spdiags([+1*e, -2*e, +1*e], [-1, 0, +1], npoints, npoints) ;

    %  2.7  Discretization matrix
    A = -nu*laplacian+c*gradient ;
    
    %  2.8  Forcing term
    forcing = M*zeros(npoints, 1) ;

    %  2.9  Dirichlet boundary conditions
    dirichlet = x == limits(1) | x == limits(2) ;
    
    Dh = distribution(dirichlet, true) ;
    Lh = distribution(dirichlet, false) ;
    
    %  2.10 System of equations
    Ah = Lh*A*Lh+Dh ;
    bh = Lh*(forcing-A*Dh*u)+Dh*u ;
    
    %  2.11 Numerical solution
    uh = Ah\bh ;
  
    
    %  3.   Data analysis

    %  3.1  Error norm
    eh = norm(M*(uh-u))/norm(M*u) ;

    %  3.2  Display
    disp(['Error norm          : ', num2str(eh)])
    disp(['Local Peclet number : ', num2str(Pel)])

    %  3.3  Analytical solution on a finer grid
    xs = linspace(0, 1, 1e+04) ;
    us = (exp(c*(xs-1)/nu)-exp(-c/nu))/(1-exp(-c/nu)) ;

    %  3.4  Plot
    figure(1) ; set(gca, 'FontName', 'Lucida Sans', 'FontSize', 14) ;
    plot(x, uh, 'bo-', xs, us, 'k-') ;
    xlabel('x') ;
    ylabel('u') ;
    title(['local Peclet number = ', num2str(Pel)]) ;
    legend('Numerical solution', 'Exact solution', 'Location', 'NorthWest') ;
    axis([0, +1, -.6, +1])
    grid on ;

end

function w = distribution(set, mode)

    if mode
        w = spdiags(set, 0, numel(set), numel(set)) ;
    else
        w = spdiags(not(set), 0, numel(set), numel(set)) ;
    end

end