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

function Unsteady_Advection_FDM_1D()

    %  1.   Parameters

    %  1.1  Domain
    length = 2e+00 ;
    
    %  1.2  Mesh
    npoints = 1e+02+1 ;

    %  1.3  Advection coefficient
    c = +1e+00 ;
    
    %  1.4  Time-integration
    to = 0e+00 ;
    tf = 2e+00 ;

    %  1.5  Theta-method
    theta = 1/2 ;
    delta_t = 1e-02 ;
    
    %  1.6  Upwind
    upwind = false ;
    
    %  1.7  Linear solver
    tolerance = 1e-06 ;
    niterations = 1e+02 ;
    solver = @(A, b, xo)bicgstabl(A, b, tolerance, niterations, [], [], xo) ;
    
    %  1.8  Eigen-value solvers
    nev = npoints-2 ;
    which = 'LM' ;
    
    %  1.9  Initial conditions
    alpha = 1e+01 ;
    sigma = 1e+00 ;
    

    %  2.   Computation

    %  2.1  Mesh points
    x = linspace(-length/2, +length/2, npoints+1)' ; 
    x = x(1:npoints) ;

    %  2.2  Mesh size
    h = length/npoints ;

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

    %  2.4  Differentiation matrices (central finite differences - 2nd order)
    derivatives{1} = (1/(2*h))*spdiags([-1*e,  0*e, +1*e], [-1, 0, +1], npoints, npoints) ;
    derivatives{2} = (1/(h^2))*spdiags([+1*e, -2*e, +1*e], [-1, 0, +1], npoints, npoints) ;
    
    %  2.5  Periodic boundary conditions
    derivatives{1}(1, end) = -1/(2*h) ;
    derivatives{1}(end, 1) = +1/(2*h) ;
    derivatives{2}(1, end) = +1/(h^2) ;
    derivatives{2}(end, 1) = +1/(h^2) ;
    
    %  2.6  Discretization matrix
    A = @(u, uo)Afun(u, uo, derivatives, c, h, upwind) ;
    
    %  2.7  Space-time discretization matrices (theta-method)
    H = @(u)Hfun(u, [], delta_t, theta, M, A) ;
    R = @(u)Rfun(u, [], delta_t, theta, M, A) ;
    
    %  2.8  Initial conditions
    uo = exp(-(alpha*x).^(2*sigma)) ;
        
    %  2.9  Numerical approximation
    ntimes = round(tf/delta_t)+1 ;
    times = (to:delta_t:tf)' ;
    
    uh = zeros(npoints, ntimes) ;
    uh(:, 1) = uo ;
    
    bar = waitbar(0, 'Time integration... please wait.') ;
    for itime = 2 : ntimes
        bh = R(uh(:, itime-1)) ;
        [uh(:, itime), flag, resvec] = solver(H, bh, uh(:, itime-1)) ;
        if flag
            warning('solver did not reach prescribed tolerance.')
        end
        waitbar(itime/ntimes, bar) ;
    end
    close(bar)
    
    %  2.10 Linear stability analysis
    Ao = Amat([], derivatives, c, h, upwind) ;
    
    [~, D] = eigs(-Ao, +M, nev, which) ;
    lambda = diag(D) ;
    
    %  2.11 Numerical stability analysis
    Ro = M/delta_t+(theta-1)*Ao ;
    Ho = M/delta_t+theta*Ao ;
    
    [~, D] = eigs(Ro, Ho, nev, which) ;
    gamma = diag(D) ;
    

    %  3.   Data analysis

    %  3.1  Error norm
    eh = zeros(ntimes, 1) ;
    u = zeros(npoints, ntimes) ;
    for itime = 1 : ntimes
        u(:, itime) = exp(-(alpha*(mod(x-c*times(itime), -length))).^(2*sigma)) ...
                    + exp(-(alpha*(mod(x-c*times(itime), +length))).^(2*sigma)) ;  
        eh(itime) = norm(M*(uh(:, itime)-u(:, itime)))/norm(M*u(:, itime)) ;
        if rem(itime-1, 10) == 0
            disp('Time step/Time/Error') ;
        end
        disp([itime, times(itime), eh(itime)]) ;
    end
        
    %  3.2  Solutions plot
    figure(1) ; clf ; 
    subplot(2, 2, 1) ;
    for itime = 2 : 2 : ntimes
        plot(x, u(:, itime), 'kx-', x, uh(:, itime), 'bo-') ;
        xlabel('x') ; ylabel('u') ;
        title(['Solutions at time t = ', num2str(times(itime))]) ;
        axis([min(x), max(x), min(uo)-0.5, max(uo)+0.5]) ;
        grid on ;
        legend('Analytical solution', 'Numerical solution')
        drawnow ;
    end
    
    %  3.3  Space-time plot of the numerical solution
    figure(1) ; subplot(2, 2, 2) ;
    
    [X, T] = ndgrid(x, times) ;
    
    pcolor(X, T, uh) ; shading interp ;
    xlabel('space') ; ylabel('time') ; 
    title('Numerical solution') ;
    caxis([min(uo), max(uo)]) ; colorbar ;
    
    %  3.4  Linear stability plot
    figure(1) ; subplot(2, 2, 3) ;
    plot(real(diag(lambda)), imag(diag(lambda)), 'bo') ;
    xlabel('Re(\lambda)') ; ylabel('Im(\lambda)') ; 
    title('Spectrum of the linear stability analysis') ;
    grid on ; axis equal ;
    
    %  3.5  Numerical stability plot
    figure(1) ; subplot(2, 2, 4) ;
    
    phi = linspace(0, 2*pi, 100) ;
    circle = exp(1i*phi) ;
    
    plot(real(diag(gamma)), imag(diag(gamma)), 'bx', real(circle), imag(circle), 'k-') ;
    xlabel('Re(\gamma)') ; ylabel('Im(\gamma)') ;
    title('Spectrum of the numerical stability analysis') ;
    grid on ; axis equal ;

end

function w = Afun(u, uo, derivatives, c, h, upwind)
    
    nuh = upwind*c*h/2 ;
    w = -nuh*derivatives{2}*u+c*derivatives{1}*u ;
    
end

function w = Rfun(u, uo, delta_t, theta, M, A)

    w = M/delta_t*u+(theta-1)*A(u, uo) ;
    
end

function w = Hfun(u, uo, delta_t, theta, M, A)

    w = M/delta_t*u+theta*A(u, uo) ;
    
end

function w = Amat(uo, derivatives, c, h, upwind)
    
    nuh = upwind*c*h/2 ;
    w = -nuh*derivatives{2}+c*derivatives{1} ;
    
end

function [x, flag, resvec] = picard(A, b, tolerance, niterations, x, solver)
  
    flag = true ;
    resvec = zeros(niterations, 1) ;
    dx = 0*x ;

    % Iterations
    normb = norm(b) ;
    for iteration = 1 : niterations
        residuals = b-A(x, x) ;
        resvec(iteration) = norm(residuals)/normb ;
        % Convergence
        if resvec(iteration) <= tolerance
            flag = false ;
            break ;
        end
        Ao = @(dx)A(dx, x) ;
        [dx, flag] = solver(Ao, residuals, 0*dx) ;
        if flag
            warning('linear solver did not reach prescribed tolerance.')
        end
        x = x+dx ;
    end
    resvec = resvec(1:iteration) ;
  
end

function u = burgers_solution(x, time, nu, nterms)

    n = 1:nterms ;
    npoints = numel(x) ;
    u = zeros(npoints, 1) ;
    
    for iterm = 1 : nterms+1
        F = @(x)exp((cos(pi*x)-1)/(2*pi*nu)).*cos((iterm-1).*pi*x);
        a(iterm) = quadgk(F, 0, 1) ;
    end
   
    for ipoint = 1 : npoints
    	u(ipoint) = ...
                (2*pi*nu*sum(exp(-(n.^2)*(pi^2)*nu*time)*2.*a(n+1).*n.*sin(n*pi*x(ipoint))))/ ...
                (a(1)+sum(exp(-(n.*n)*(pi^2)*nu*time)*2.*a(n+1).*cos(n*pi*x(ipoint))));           
    end
    
end