% -------------------------------------------------------------------------------------- %
% Methodes de discretisation en fluides                                                  %
% M.A. Habisreutinger                                                                    %
% Jeudi 9 avril 2019                                                                     %
% -------------------------------------------------------------------------------------- %
%                                                                                        %
%                            SOLUTION IV                                                 %
%                                                                                        %
% -------------------------------------------------------------------------------------- %

function SteadyPoiseuile2D_Piso()


    %% 1.   PARAMETERS

    %  1.1  Mesh
    ndimensions = 2 ;
    npoints = [32, 16] ;
    L = [8, 1] ;

    %  1.2  Fluid density and viscosity
    rho = 1e+00; 
    mu = 1e+00 ;
    
    %  1.3  Body force
    g.x = 0e+00 ;
    g.y = 0e+00 ;

    %  1.4  Pressure drop
    dpdx = -1e+00 ;

    %  1.5  Velocity-profile plot
    location = 1/4 ;
    lambda = 5e+00 ;

    %  1.6  PISO parameters
    niterations = 2e+00 ;
    tolerances.momentum = 1e-06 ;
    tolerances.mass = 1e-06 ;
    

    %% 2.   COMPUTATION

    %  2.1  Number of nodes
    nnodes = prod(npoints) ;

    %  2.2  Mesh coordinates and size
    for idim = 1 : ndimensions
      xi{idim} = linspace(-L(idim)/2, +L(idim)/2, npoints(idim)) ;
      h(idim) = L(idim)/(npoints(idim)-1) ;
    end
    
    [X, Y] = ndgrid(xi{1}, xi{2}) ;
    coordinates.x = reshape(X, nnodes, 1) ;
    coordinates.y = reshape(Y, nnodes, 1) ;

    %  2.3  Mesh boundaries
    mask = sparse(npoints(1), npoints(2)) ;
    mask(1, :) = true ;
    boundaries.left = reshape(mask, nnodes, 1) ;

    mask = sparse(npoints(1), npoints(2)) ;
    mask(end, :) = true ;
    boundaries.right = reshape(mask, nnodes, 1) ;

    mask = sparse(npoints(1), npoints(2)) ;
    mask(:, 1) = true ;
    boundaries.bottom = reshape(mask, nnodes, 1) ;

    mask = sparse(npoints(1), npoints(2)) ;
    mask(:, end) = true ;
    boundaries.top = reshape(mask, nnodes, 1) ;
    
    %  2.4  Mesh normals at boundaries
    normals.left.x = -boundaries.left ;
    normals.left.y = zeros(nnodes, 1) ;
    
    normals.right.x = +boundaries.right ;
    normals.right.y = zeros(nnodes, 1) ;
    
    normals.bottom.x = zeros(nnodes, 1) ;
    normals.bottom.y = -boundaries.bottom ;
    
    normals.top.x = zeros(nnodes, 1) ;
    normals.top.y = +boundaries.top ;
    
    %  2.5  Discrete differential operators (collocated fintes differences)
    for idim = 1 : ndimensions
      % 2.5.1 Identity
      e = ones(npoints(idim), 1) ;
      I{idim} = spdiags(e, 0, npoints(idim), npoints(idim)) ;
      % 2.5.2 Laplacian with Neumann boundary conditions
      Ai{idim} = spdiags((1/h(idim)^2)*[+1*e -2*e +1*e], [-1, 0, +1], npoints(idim), npoints(idim)) ;
      Ai{idim}(1, 1) = +1/h(idim) ; Ai{idim}(1, 2) = -1/h(idim) ;
      Ai{idim}(end, end) = +1/h(idim) ; Ai{idim}(end, end-1) = -1/h(idim) ;
      % 2.5.3 Gradient (2nd order and shifted scheme at boundaries)
      Gi{idim} = spdiags((1/(2*h(idim)))*[-1*e  0*e +1*e], [-1, 0, +1], npoints(idim), npoints(idim)) ;
      Gi{idim}(1, 1) = -1/h(idim) ;
      Gi{idim}(1, 2) = +1/h(idim) ;
      Gi{idim}(end, end-1) = -1/h(idim) ;
      Gi{idim}(end, end) = +1/h(idim) ;
      % 2.5.4 Laplacian (2nd order and shifted scheme at boundaries)
      Li{idim} = spdiags((1/h(idim)^2)*[+1*e -2*e +1*e], [-1, 0, +1], npoints(idim), npoints(idim)) ;
      Li{idim}(1, 1) = +1/h(idim)^2 ;
      Li{idim}(1, 2) = -2/h(idim)^2 ;
      Li{idim}(1, 3) = +1/h(idim)^2 ;
      Li{idim}(end, end-2) = +1/h(idim)^2 ;
      Li{idim}(end, end-1) = -2/h(idim)^2 ;
      Li{idim}(end, end) = +1/h(idim)^2 ;
    end

    %  2.6  Tensorisation for a 2D problem
    Ao = kron(I{2}, Ai{1})+kron(Ai{2}, I{1}) ;
    laplacian = kron(I{2}, Li{1})+kron(Li{2}, I{1}) ;
    gradient.x = kron(I{2}, Gi{1}) ;
    gradient.y = kron(Gi{2}, I{1}) ;

    %  2.7  Source term
    source.x = g.x*ones(nnodes, 1) ;
    source.y = g.y*ones(nnodes, 1) ;
    
    %  2.8  Exact solution
    velocity.x = (1/2)*dpdx*(coordinates.y-L(2)/2).*(coordinates.y+L(2)/2) ;
    velocity.y = zeros(nnodes, 1) ;
    pressure = dpdx*coordinates.x ;

    %  2.9  Velocity guess for PISO
    vh.x = zeros(nnodes, 1) ;
    vh.y = zeros(nnodes, 1) ;
    ph = zeros(nnodes, 1) ;
    
    %  2.10 PISO boundary conditions
    
    %  2.10.1 Velocity x-component
    dirichlet.x = boundaries.bottom | boundaries.top ;
    
    %  2.10.2 Velocity y-component
    dirichlet.y = boundaries.bottom | boundaries.top | boundaries.left | boundaries.right ;
    
    %  2.10.3 Velocity pressure
    dirichlet.p = boundaries.left | boundaries.right ;
    neumann.p = boundaries.bottom | boundaries.top ;
    
    %  2.11 PISO operators
    
    %  2.11.1 Velocity x-component
    D = spdiags(dirichlet.x, 0, nnodes, nnodes) ;
    F = spdiags(~dirichlet.x, 0, nnodes, nnodes) ;
    Piso.x = F*(-mu)*Ao+D ;
    
    %  2.11.2 Velocity y-component
    D = spdiags(dirichlet.y, 0, nnodes, nnodes) ;
    F = spdiags(~dirichlet.y, 0, nnodes, nnodes) ;
    Piso.y = F*(-mu)*Ao+D ;
    
    %  2.11.3 Pressure
    D = spdiags(dirichlet.p, 0, nnodes, nnodes) ;
    F = spdiags(~dirichlet.p, 0, nnodes, nnodes) ;
    Piso.p = F*Ao+D ;
    
    %  2.12 PISO iterations
    for iteration = 1 : niterations
        
        %  2.12.1 Residuals
        residuals.momentum_x(iteration) = norm(-mu*laplacian*velocity.x+gradient.x*ph+rho*source.x) ;
        residuals.momentum_y(iteration) = norm(-mu*laplacian*velocity.y+gradient.y*ph+rho*source.y) ;
        residuals.mass(iteration) = norm(gradient.x*velocity.x+gradient.y*velocity.y) ;
        
        %  2.12.2 Convergence
        booleans(1) = residuals.momentum_x(iteration) < tolerances.momentum ;
        booleans(2) = residuals.momentum_y(iteration) < tolerances.momentum ;
        booleans(3) = residuals.mass(iteration) < tolerances.mass ;
        
        if all(booleans)
            disp(['√ PISO converged to prescribed tolerances afer ', num2str(iteration-1),' iteration(s).'])
            break ;
        elseif iteration == niterations
            disp(['x PISO did not converge to prescribed tolerances afer ', num2str(niterations),' iteration(s)!']) ;
            break ;
        end
        
        %  2.12.3 Pressure update
        grad_ph = grad_p(rho, mu, vh, source, laplacian) ;
        dph_dn = grad_ph.x.*normals.bottom.x+grad_ph.y.*normals.bottom.y ...
               + grad_ph.x.*normals.top.x+grad_ph.y.*normals.top.y ; 

        divg = rho*(gradient.x*source.x+gradient.y*source.y) ;
        rhs = not(dirichlet.p).*divg+dirichlet.p.*pressure+neumann.p.*dph_dn ;
        ph = Piso.p\rhs ;

        %  2.12.4 Velocity x-component update
        rhs = not(dirichlet.x).*(-gradient.x*ph+rho*source.x)+dirichlet.x.*velocity.x ;
        vh.x = Piso.x\rhs ;

        %  2.12.5 Velocity y-component update
        rhs = not(dirichlet.y).*(-gradient.y*ph+rho*source.y)+dirichlet.y.*velocity.y ;
        vh.y = Piso.y\rhs ;
    end
    

    %% 3.   DATA ANALYSIS
    
    %  3.1  Error norms
    errors.velocity_x = norm(vh.x-velocity.x) ;
    errors.velocity_y = norm(vh.y-velocity.y) ;
    errors.pressure = norm(ph-pressure) ;
    errors
    
    %  3.2  Piso residuals
    figure(1) ;
    semilogy(1:iteration, residuals.momentum_x+eps, 'LineWidth', 2) ; hold on ;
    semilogy(1:iteration, residuals.momentum_y+eps, 'LineWidth', 2) ;
    semilogy(1:iteration, residuals.mass+eps, 'LineWidth', 2) ; hold off ;
    grid on ;
    xlabel('iterations') ; ylabel('residuals') ; title('Piso residuals') ;
    legend('momentum-x', 'momentum-y', 'mass')
    
    %  3.3  Solution plots
    Vxh = reshape(vh.x, npoints) ; 
    Vyh = reshape(vh.y, npoints) ; 
    ix = floor(npoints(1)*location) ;

    %  3.3.1 x-velocity plot
    figure(2) ; subplot(3, 1, 1) ; 
    pcolor(X, Y, Vxh) ; hold on ;
    quiver(X(ix, :), Y(ix, :), Vxh(ix, :), Vyh(ix, :), lambda) ; hold off ;
    shading interp ; colorbar ; colormap('hot') ; axis equal ; axis tight ;
    xlabel('x') ; ylabel('y') ; title('Approximated x-velocity field') ;

    %  3.3.2 y-velocity plot
    figure(2) ; subplot(3, 1, 2) ; 
    pcolor(X, Y, Vyh) ; 
    shading interp ; colorbar ; colormap('hot') ; axis equal ; axis tight ;
    xlabel('x') ; ylabel('y') ; title('Approximated y-velocity field') ;

    %  3.3.3 pressure error plot
    figure(2) ; subplot(3, 1, 3) ; 
    pcolor(X, Y, reshape(ph-pressure, npoints)) ; 
    shading interp ; colorbar ; colormap('hot') ; axis equal ; axis tight ;
    xlabel('x') ; ylabel('y') ; title('Error on pressure field') ;

end

% -------------------------------------------------------------------------------------- %

function w = grad_p(rho, mu, velocity, source, laplacian)

    w.x = mu*laplacian*velocity.x+rho*source.x ;
    w.y = mu*laplacian*velocity.y+rho*source.y ;

end

% -------------------------------------------------------------------------------------- %