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

function Sol01()

    % Convergence study
    nrefinements = 8 ;
    for ir = 1 : nrefinements
        npoints(ir)  = 2^(ir+2)+1 ;
        [eh(ir, 1), h(ir, 1)] = ErrorDiffusion1D(2, npoints(ir)) ;
        [eh(ir, 2), ~] = ErrorDiffusion1D(4, npoints(ir)) ;
    end

    % Convergence plot
    figure(1) ;
    set(gca, 'FontName', 'Lucida Sans', 'FontSize', 14) ;
    loglog(h, eh(:, 1), 'r', h, (3*h).^2, 'r:', h, eh(:, 2), 'b', h, (2.4*h).^4, 'b:') ;
    xlabel('mesh size') ; ylabel('L2 norm of error') ;
    grid on ;
    title('Convergence study') ;
    legend('2^{nd}-order scheme', 'O(h^2)', '4^{th}-order scheme', 'O(h^4)', 'Location', 'Best') ;

end

function [eh, h] = ErrorDiffusion1D(order, npoints)

    %  1.   Parameters

    %  1.1  Domain
    limits = [-1, +1] ;

    %  1.2  Diffusion coefficient
    nu = 1e+00 ;


    %  2.   Computation

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

    %  2.2  Mesh normals
    normals = zeros(npoints, 1) ;
    normals(1) = -1 ;
    normals(end) = +1 ;
    
    %  2.3  Mesh size
    h = (limits(2)-limits(1))/(npoints-1) ;

    %  2.4  Analytical solution
    k = 2*pi/(limits(2)-limits(1)) ;
    u = cos(k*x+pi/4) ;
    u_x = -k*sin(k*x+pi/4) ;

    %  2.5  Mass matrix (colocation method)
    e = ones(npoints, 1) ;
    M = spdiags(e, 0, npoints, npoints) ;
    
    %  Correction for boundary conditions
    M(1, 1) = 0 ;
    M(end, end) = 0 ;

    %  2.6  Laplacian matrix
    if order == 2
        laplacian = DifferentiationMatrix(2, -1:+1, npoints, h) ;
        %  = (1/h^2)*spdiags([+1*e, -2*e, +1*e], -1:+1, npoints, npoints) ;
        
        %  Correction for Neumann boundary conditions
        taylor = TaylorCoefficients(0:+2, h) ;
        laplacian(1, 1) = +taylor(2, 1) ; % = +3/(2*h)
        laplacian(1, 2) = +taylor(2, 2) ; % = -4/(2*h)
        laplacian(1, 3) = +taylor(2, 3) ; % = +1/(2*h)

        taylor = TaylorCoefficients(-2:0, h) ;
        laplacian(end, end-2) = -taylor(2, 1) ; % = -1/(2*h)
        laplacian(end, end-1) = -taylor(2, 2) ; % = +4/(2*h)
        laplacian(end, end)   = -taylor(2, 3) ; % = -3/(2*h)
    elseif order == 4
        laplacian = DifferentiationMatrix(2, -2:+2, npoints, h) ;
        % = (1/h^2)*spdiags([-1*e/12, +4*e/3, -5*e/2, +4*e/3, -1*e/12], -2:+2, npoints, npoints) ;
        
        % Shifted schemes for nodes 2 and npoints-1
        taylor = TaylorCoefficients(-1:4, h) ;
        laplacian(2, 1) = +taylor(3, 1) ; % = (1/h^2)*(+5/6.)
        laplacian(2, 2) = +taylor(3, 2) ; % = (1/h^2)*(-5/4.)
        laplacian(2, 3) = +taylor(3, 3) ; % = (1/h^2)*(-1/3.)
        laplacian(2, 4) = +taylor(3, 4) ; % = (1/h^2)*(+7/6.)
        laplacian(2, 5) = +taylor(3, 5) ; % = (1/h^2)*(-1/2.)
        laplacian(2, 6) = +taylor(3, 6) ; % = (1/h^2)*(+1/12)
        
        taylor = TaylorCoefficients(-4:1, h) ;
        laplacian(end-1, end-5) = +taylor(3, 1) ; % = (1/h^2)*(+1/12) ;
        laplacian(end-1, end-4) = +taylor(3, 2) ; % = (1/h^2)*(-1/2.) ;
        laplacian(end-1, end-3) = +taylor(3, 3) ; % = (1/h^2)*(+7/6.) ;
        laplacian(end-1, end-2) = +taylor(3, 4) ; % = (1/h^2)*(-1/3.) ;
        laplacian(end-1, end-1) = +taylor(3, 5) ; % = (1/h^2)*(-5/4.) ;
        laplacian(end-1, end)   = +taylor(3, 6) ; % = (1/h^2)*(+5/6.) ;
        
        %  Correction for Neumann boundary conditions
        taylor = TaylorCoefficients(0:4, h) ;
        laplacian(1, 1) = +taylor(2, 1) ; % = -25/(12*h) ; 
        laplacian(1, 2) = +taylor(2, 2) ; % = +48/(12*h) ;
        laplacian(1, 3) = +taylor(2, 3) ; % = -36/(12*h) ;
        laplacian(1, 4) = +taylor(2, 4) ; % = +16/(12*h) ;
        laplacian(1, 5) = +taylor(2, 5) ; % = -3./(12*h) ;
        
        taylor = TaylorCoefficients(-4:0, h) ;
        laplacian(end, end-4) = -taylor(2, 1) ; % = -3./(12*h) ;
        laplacian(end, end-3) = -taylor(2, 2) ; % = +16/(12*h) ;
        laplacian(end, end-2) = -taylor(2, 3) ; % = -36/(12*h) ;
        laplacian(end, end-1) = -taylor(2, 4) ; % = +48/(12*h) ;
        laplacian(end, end)   = -taylor(2, 5) ; % = -25/(12*h) ;
    end
    
    %  2.7  Discretization matrix
    A = -nu*laplacian ;
    
    %  2.8  Forcing term
    forcing = M*norm(k)^2*nu*u ;
    
    %  2.9  Neumann boundary conditions
    neumann = x == limits(2) ;
    forcing = forcing+neumann.*(nu*u_x.*normals) ;
    
    %  2.10 Dirichlet boundary conditions
    dirichlet = x == limits(1) ;
    
    Dh = distribution(dirichlet, true) ;
    Lh = distribution(dirichlet, false) ;
    
    %  2.11 System of equations
    Ah = Lh*A*Lh+Dh ;
    bh = Lh*(forcing-A*Dh*u)+Dh*u ;

    %  2.12 Numerical solution
    uh = Ah\bh ;
    

    %  3.   Data analysis

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

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

function w = TaylorCoefficients(stencil, h)
    
    nstencil = numel(stencil) ;
    T = zeros(nstencil, nstencil) ;
    
    for ip = 0 : nstencil-1
        for jp = 0 : nstencil-1 
            T(ip+1, jp+1) = (stencil(ip+1)*h)^jp/factorial(jp) ;
        end
    end
    
    w = inv(T) ;
    
end

function w = DifferentiationMatrix(order, stencil, npoints, h)

    e = ones(npoints, 1) ;
    taylor = TaylorCoefficients(stencil, h) ;
    
    diagonals = zeros(npoints, numel(stencil)) ;
    for idiag = 1 : numel(stencil)
        diagonals(:, idiag) = e*taylor(order+1, idiag) ;
    end
    
    w = spdiags(diagonals, stencil, npoints, npoints) ;
    
end