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

function Unsteady_Burgers_FDM_1D()

    %  1.   Parameters

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

    %  1.3  Diffusion coefficient
    nu = 1e-02 ;
    
    %  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 = 100 ;
    
    %  1.9  Non-linear solver
    nl_solver = @(A, b, x)picard(A, b, tolerance, niterations, x, solver) ;
    
    %  1.10 Analytical solution
    nterms = 5e+01 ;
    

    %  2.   Computation

    %  2.1  Mesh points
    x = linspace(0, length, 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, nu, h, upwind) ;
    
    %  2.7  Space-time discretization matrices (theta-method)
    H = @(u, uo)Hfun(u, uo, delta_t, theta, M, A) ;
    R = @(u, uo)Rfun(u, uo, delta_t, theta, M, A) ;
    
    %  2.8  Initial conditions
    uo = sin(pi*x) ;

    %  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-1)) ;
        [uh(:, itime), flag, resvec] = nl_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 at initial conditions
    Lo = Lmat(uo, derivatives, nu, h, upwind) ;
    
    [~, D] = eigs(-Lo, +M, nev, 'LR') ;
    lambda = diag(D) ;
        
    %  2.11 Numerical stability analysis at initial conditions
    Ro = M/delta_t+(theta-1)*Lo ;
    Ho = M/delta_t+theta*Lo ;
    
    [~, D] = eigs(Ro, Ho, nev, 'LM') ;
    gamma = diag(D) ;
        
        
    %  3.   Data analysis

    %  3.1  Error norm
    eh = zeros(ntimes, 1) ;
    u = zeros(npoints, ntimes) ;
    bar = waitbar(0, 'Computing the exact solution... please wait.') ;
    for itime = 1 : ntimes
        u(:, itime) = burgers_solution(x, times(itime), nu, nterms) ;
        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)]) ;
        waitbar(itime/ntimes, bar) ;
    end
    close(bar)
    
    %  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)), 'bo', 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, nu, h, upwind)
    
    nuh = upwind*max(abs(uo))*h/2 ;
    w = -(nu+nuh)*(derivatives{2}*u)+derivatives{1}*(uo.*u/2) ;
    
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 = Lmat(uo, derivatives, nu, h, upwind)
    
    nuh = upwind*max(abs(uo))*h/2 ;
    
    Uo = spdiags(uo, 0, numel(uo), numel(uo)) ;
    Uox = spdiags(derivatives{1}*uo, 0, numel(uo), numel(uo)) ;
    
    w = -(nu+nuh)*derivatives{2}+Uo*derivatives{1}+Uox ;
    
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