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

function Unsteady_Advection_FDM_2D()

    %  1.   Parameters

    %  1.1  Domain
    ndimensions = 2 ;
    lengths = [1e+00, 1e+00] ;
    
    %  1.2  Mesh
    npoints = 2^7*lengths ; 

    %  1.3  Advection coefficient
    c = [+1e+00, +1e+00] ;
    
    %  1.4  Time-integration
    to = 0e+00 ;
    tf = 1e+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  Initial conditions
    alpha = 1e+01 ;
    sigma = 1e+00 ;
    

    %  2.   Computation

    %  2.1  Nodal distributions
    for idim = 1 : ndimensions
        xi{idim} = linspace(-lengths(idim)/2, +lengths(idim)/2, npoints(idim)+1)' ; 
        xi{idim} = xi{idim}(1:npoints(idim)) ;
        h(idim) = lengths(idim)/npoints(idim) ;
    end

    %  2.2  Mesh points
    [X{1}, X{2}] = ndgrid(xi{1}, xi{2}) ;
    for idim = 1 : ndimensions
        x(:, idim) = reshape(X{idim}, prod(npoints),1) ;
    end

     %  2.4  Discrete operators
    for idim = 1 : ndimensions
        e = ones(npoints(idim), 1) ;
        % Mass matrix (colocation method)
        Mi{idim} = spdiags(e, 0, npoints(idim), npoints(idim)) ;
        % Differentiation matrices (central finite differences - 2nd order)
        Di{1}{idim} = (1/(2*h(idim)))*spdiags([-1*e, 0*e, +1*e], [-1, 0, +1], npoints(idim), npoints(idim)) ;
        Di{1}{idim}(1, end) = -1/(2*h(idim)) ;
        Di{1}{idim}(end, 1) = +1/(2*h(idim)) ;
        Di{2}{idim} = 1/h(idim)^2*spdiags([+1*e, -2*e, +1*e], [-1, 0, +1], npoints(idim), npoints(idim)) ;
        Di{2}{idim}(1, end) = +1/(h(idim)^2) ;
        Di{2}{idim}(end, 1) = +1/(h(idim)^2) ;
    end
    M = kron(Mi{2}, Mi{1}) ;
    
    gradient{1} = kron(Mi{2}, Di{1}{1}) ;
    gradient{2} = kron(Di{1}{2}, Mi{1}) ;
    
    laplacian = kron(Mi{2}, Di{2}{1})+kron(Di{2}{2}, Mi{1}) ;
    
    %  2.6  Discretization matrix
    A = @(u, uo)Afun(u, uo, gradient, laplacian, 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.7 Initial conditions
    [~, r] = cart2pol(x(:, 1), x(:, 2)) ;
    uo = exp(-(alpha*r).^(2*sigma)) ;
        
    %  2.8  Numerical approximation
    ntimes = round(tf/delta_t)+1 ;
    times = (to:delta_t:tf)' ;

    uh = zeros(prod(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)
    

    %  3.   Data analysis

    %  3.2  Numerical solution plot
    figure(1) ;
    for itime = 1 : 2 : ntimes
        pcolor(reshape(x(:, 1), npoints), reshape(x(:, 2), npoints), reshape(uh(:, itime), npoints)) ; 
        shading interp ; colorbar ;
        caxis([0, 1]) ; colormap('hot') ;
        axis([-lengths(1)/2, +lengths(1)/2, -lengths(2)/2, +lengths(2)/2]) ;
        xlabel('x') ; ylabel('y') ; title('Numerical solution') ;
        drawnow ;
    end

end

function w = Afun(u, uo, gradient, laplacian, c, h, upwind)
    
    nuh = upwind*norm(c)*norm(h)/2 ;
    w = -nuh*(laplacian*u)+c(1)*(gradient{1}*u)+c(2)*(gradient{2}*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