% -------------------------------------------------------------------------------------- %
% Methodes de discretisation en fluides                                                  %
% Laboratoire II - Solution                                                              %
% M.A. Habisreutinger                                                                    %
% Jeudi 28 mars 2024                                                                     %
% -------------------------------------------------------------------------------------- %

function Sol02()
  
    %  1.  Mesh
    npoints = 16 ;

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

    %  3.  Refinement levels
    nrefinements = 16 ;

    %  4.  Comparing Centered and Upwind shemes
    figure(1) ;
    Comparison(npoints, nu, c) ;

    %  5.  Convergence study
    figure(2) ;
    Convergence(nrefinements, nu, c) ;

end

function Comparison(npoints, nu, c)
  
    %  1.  Analytical solution on a fine grid
    x = linspace(0, +1, 1e+03) ;
    u = (exp(c*(x-1)/nu)-exp(-c/nu))/(1-exp(-c/nu)) ;

    %  2.  Numerical solutions (centered and upwind schemes)
    [~ , uh, ~, ~  ] = AdvectionDiffusion1D(npoints, nu, c, 'centered') ;
    [~ , us, ~, ~  ] = AdvectionDiffusion1D(npoints, nu, c, 'upwind_2nd_order') ;
    [xh, up, ~, Pel] = AdvectionDiffusion1D(npoints, nu, c, 'upwind_1st_order') ;

    %  3.  Plot
    plot(x, u, 'k-', xh, uh, 'mo-', xh, us, 'cd-', xh, up, 'kx-') ;
    xlabel('x') ; ylabel('u') ;
    title(['Solutions of the problem,  local Peclet number = ', num2str(Pel)]) ;
    legend('Exact solution', 'Centered scheme', 'Upwind 2nd-order',  ...
        'Upwind 1st-order', 'Location', 'NorthWest') ;
    grid on ;
  
end

function Convergence(nrefinements, nu, c)
  
    %  1.  Error analysis
    npoints = zeros(nrefinements, 1) ;
    eh = zeros(nrefinements, 1) ;
    es = zeros(nrefinements, 1) ;
    ep = zeros(nrefinements, 1) ;
    Pel = zeros(nrefinements, 1) ;

    bar = waitbar(0, 'Studying convergence... please wait.') ;
    for ir = 1 : nrefinements
        npoints(ir) = 2^(ir+1)+1 ;
        [~, ~, eh(ir), ~      ] = AdvectionDiffusion1D(npoints(ir), nu, c, 'centered') ;
        [~, ~, es(ir), ~      ] = AdvectionDiffusion1D(npoints(ir), nu, c, 'upwind_2nd_order') ;
        [~, ~, ep(ir), Pel(ir)] = AdvectionDiffusion1D(npoints(ir), nu, c, 'upwind_1st_order') ;
        waitbar(ir/nrefinements, bar) ;
    end
    close(bar) ;

    %  2.  Plot
    loglog(Pel, eh, 'mo:', Pel, .15*Pel.^2, 'm-',  ...
        Pel, es, 'cd:', Pel, .80*Pel.^2, 'c-',  ...
        Pel, ep, 'kx:', Pel, .50*Pel.^1, 'k-') ;
    xlabel('Local Peclet number') ;
    ylabel('Error') ;
    legend('Centered scheme', 'O(h^2)', 'Upwind 2nd-order', 'O(h^2)',  ...
        'Upwind 1st-order', 'O(h^1)', 'Location', 'Best') ;
    grid on ;

end

function [x, uh, eh, Pel] = AdvectionDiffusion1D(npoints, nu, c, method)

    %  1.   Parameters

    %  1.1  Mesh
    limits = [0, +1] ;
    
    
    %  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
    if strcmp(method, 'upwind_2nd_order')
        gradient = (1/(2*h))*spdiags([+1*e, -4*e, +3*e], [-2, -1, 0], npoints, npoints) ;
        gradient(2, :) = 0 ;
        gradient(2, 1) = -1/h ;
        gradient(2, 2) = +1/h ;
    else
        gradient = (1/(2*h))*spdiags([-1*e,  0*e, +1*e], [-1, 0, +1], npoints, npoints) ;
    end
    laplacian = (1/(h^2))*spdiags([+1*e, -2*e, +1*e], [-1, 0, +1], npoints, npoints) ;

    %  2.7  Discretization matrix
    nuh = nu*artificial_diffusion(Pel, method) ;
    A = -(nu+nuh)*laplacian+c*gradient ;

    %  2.8  Forcing term
    forcing = M*zeros(npoints, 1) ;

    %  2.9  Dirichlet boundary conditions
    dirichlet = x == 0 | x == 1 ;

    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 ;

    %  2.12 Error norms
    eh = norm(M*(uh-u))/norm(M*u) ;
  
end

function w = artificial_diffusion(Pel, method)
  
    switch method
      case 'centered'
        w = 0 ;
      case 'upwind_1st_order'
        w = Pel ;
      case 'sharfeter_gummel'
        w = Pel-1+(2*Pel)/(exp(2*Pel)-1) ;
      case 'upwind_2nd_order'
        w = 0 ;
    end

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