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

function Unsteady_Advection_FDM_3D()

    %  1.   Parameters

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

    %  1.3  Advection coefficient
    c = [+0e+00, +1e+00, +0e+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 ;
    
    %  1.9  Vtk output
    output.vtk = false ;
    output.filename = 'Unsteady_Advection_FDM_3D' ;
    

    %  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}, X{3}] = ndgrid(xi{1}, xi{2}, xi{3}) ;
    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{3}, kron(Mi{2}, Mi{1})) ;
    
    gradient{1} = kron(Mi{3}, kron(Mi{2}, Di{1}{1})) ;
    gradient{2} = kron(Mi{3}, kron(Di{1}{2}, Mi{1})) ;
    gradient{3} = kron(Di{1}{3}, kron(Mi{2}, Mi{1})) ;
    
    laplacian = kron(Mi{3}, kron(Mi{2}, Di{2}{1})) ...
              + kron(Mi{3}, kron(Di{2}{2}, Mi{1})) ...
              + kron(Di{2}{3}, kron(Mi{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] = cart2sph(x(:, 1), x(:, 2), x(:, 3)) ;
    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.1  Vtk output
    if output.vtk 
        bar = waitbar(0, 'Writing vtk output... please wait.') ;
        for itime = 1 : ntimes
            filename = strcat(output.filename, '_', num2str(itime)) ;
            WriteVtkOutput(filename, xi, npoints, uh(:, itime)) ;
            waitbar(itime/ntimes, bar) ;
        end
        close(bar)
    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

function WriteVtkOutput(filename, xi, npoints, uh)

    % Open file
    file = fopen(strcat(filename, '.vtk'), 'w') ;

    % Header
    fprintf(file, '%s\n', '# vtk DataFile Version 3.1') ;
    fprintf(file, '%s\n', 'Export Q') ;
    fprintf(file, '%s\n', 'ASCII') ;
    fprintf(file, '%s\n', 'DATASET RECTILINEAR_GRID') ;

    % Mesh
    fprintf(file, '%s ', 'DIMENSIONS') ;
    fprintf(file, '%u %u %u\n', npoints) ;

    fprintf(file, '%s ', 'X_COORDINATES') ;
    fprintf(file, '%u ', npoints(1)) ;
    fprintf(file, '%s\n', 'FLOAT') ;

    fprintf(file, '%15.5f', xi{1}) ; fprintf(file, '\n') ;

    fprintf(file, '%s ', 'Y_COORDINATES') ;
    fprintf(file, '%u ', npoints(2)) ;
    fprintf(file, '%s\n', 'FLOAT') ;

    fprintf(file, '%15.5f', xi{2}) ; fprintf(file,'\n') ;

    fprintf(file, '%s ', 'Z_COORDINATES') ;
    fprintf(file, '%u ', npoints(3)) ;
    fprintf(file, '%s\n', 'FLOAT') ;

    fprintf(file, '%15.5f', xi{3}) ; fprintf(file,'\n') ;

    fprintf(file, '%s ', 'POINT_DATA') ;
    fprintf(file, '%u\n', prod(npoints)) ;

    % Fields
    fprintf(file, '%s ', 'SCALARS') ;
    fprintf(file, '%s ', 'uh') ;
    fprintf(file, '%s\n', 'FLOAT') ;
    fprintf(file, '%s\n', 'LOOKUP_TABLE default') ;
    fprintf(file, '%15.5f', uh') ; fprintf(file, '\n') ;

    % Close file
    fclose(file) ;
  
end