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

clear all ; close all ; 
format compact ; format short e ;


%% 1.   PARAMETERS

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

%  1.2  Fluid viscosity
mu = 1e+00 ;

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

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


%% 2.   COMPUTATION

%  2.1  Number of variables, nodes and unknowns
nvars = ndimensions+1 ;
nnodes = prod(npoints) ;
nuks = nvars*nnodes ;

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

%  2.3  Mesh
[X, Y] = ndgrid(xi{1}, xi{2}) ;
x = reshape(X, nnodes, 1) ;
y = reshape(Y, nnodes, 1) ;

%  2.4  Boundary descriptors
mask = sparse(npoints(1), npoints(2)) ;
mask(1, :) = true ;
left = reshape(mask, nnodes, 1) ;

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

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

mask = sparse(npoints(1), npoints(2)) ;
mask(:, end) = true ;
top = reshape(mask, nnodes, 1) ;

%  2.5  Discrete operators
for idim = 1 : ndimensions
  e = ones(npoints(idim), 1) ;
  % Central finite differences (2nd order)
  Mi{idim} = spdiags(e, 0, npoints(idim), npoints(idim)) ;
  Ai{idim} = spdiags(-mu*(1/h(idim)^2)*[+1*e -2*e +1*e], [-1, 0, +1], npoints(idim), npoints(idim)) ;
  Di{idim} = spdiags(+ (1/(2*h(idim)))*[-1*e  0*e +1*e], [-1, 0, +1], npoints(idim), npoints(idim)) ;
  % Neumann boundary conditions
  Ai{idim}(1, 1) = -mu*(-1/h(idim)^2) ;
  Ai{idim}(end, end) = -mu*(-1/h(idim)^2) ;
  % Shifted scheme (1st order)
  Di{idim}(1, 1) = -1/h(idim) ;
  Di{idim}(1, 2) = +1/h(idim) ;
  Di{idim}(end, end-1) = -1/h(idim) ;
  Di{idim}(end, end) = +1/h(idim) ;
end

%  2.6  Tensorisation for a 2D problem
Mo = kron(Mi{2}, Mi{1}) ;
Ao = kron(Mi{2}, Ai{1})+kron(Ai{2}, Mi{1}) ;
Dx = kron(Mi{2}, Di{1}) ;
Dy = kron(Di{2}, Mi{1}) ;
O = sparse(nnodes, nnodes) ;

%  2.7  Assembling global matrices
M = [Mo, O,  O; ...
      O, Mo, O; ...
      O, O,  O] ;

A = [Ao, O,  Dx; ...
      O, Ao, Dy; ...
     Dx, Dy, O] ;

%  2.8  Forcing term
f = zeros(nuks, 1) ;
Mf = M*f ;

%  2.9  Exact solution
vx = (1/2)*dpdx*(y-L(2)/2).*(y+L(2)/2) ;
vy = zeros(nnodes, 1) ;
p = dpdx*x ;

u = [vx; vy; p] ;

%  2.10 Correction for Dirichlet boundary conditions
c = [(bottom | top); ...
  +  (bottom | top | left | right); ...
  +  (left | right)] ;

C = spdiags( c, 0, nuks, nuks) ;
F = spdiags(~c, 0, nuks, nuks) ;

Mf = F*Mf+C*u ;
A = F*A+C ;

%  2.13 Solution of the discrete equations (direct method)
uh = A\Mf ;

%  2.14 Retrieving the variables
indices = 1:nnodes ;
vxh = uh(indices) ;
vyh = uh(indices+nnodes) ;
ph = uh(indices+2*nnodes) ;


%% 3.   DATA ANALYSIS

%  3.1  Error norms
eh = norm(M*(uh-u))/norm(M*u) ;

%  3.2  Display
disp(['Error norm : ', num2str(eh)]) ;

%  3.3  Global matrices plot
figure(1) ;
spy(A) ;
xlabel('j') ; ylabel('i') ; 
title('Structure of the global discretization matrix') ;

%  3.4  Solution plots
Vxh = reshape(vxh, npoints) ; 
Vyh = reshape(vyh, npoints) ; 
ix = floor(npoints(1)*location) ;

%  3.4.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.4.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.4.3 pressure error plot
figure(2) ; subplot(3, 1, 3) ; 
pcolor(X, Y, reshape(ph-p, npoints)) ; 
shading interp ; colorbar ; colormap('hot') ; axis equal ; axis tight ;
xlabel('x') ; ylabel('y') ; title('Error on pressure field') ;

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