% ------------------------------------------------------------------------
% 
% Problem: 2D steady-state advection-diffusion, harmonic RHS:
%         div(rho*w*phi) = div(Gamma*grad(phi)) + sin(x) * sin(y)
%
% Domain: square [0,2pi] x [0,2pi]
%
% Grid: structured, cartesian, constant spacing, dx = dy
%
% Boundary conditions (BC):
%   Dirichet ( phi           = phi_BC ) on the lower side ( y = 0 )
%   Neumann  ( Gamma * phi_n = Jd_BC  ) on the other sides
%
% Discretization method: Finite Volume (FV)
%
% Solution methodology: direct and Conjugate Gradient 
%                       (sparse coefficient matrix)
%
% Date: 30/3/2016
% Author: R. Zamolo
%
% ------------------------------------------------------------------------
%
% Domain, coordinates systems and boundary conditions:
%
%         ^ y( j ) 
%         |
%         |
%         |    Neumann
%     2pi +---------------+
%         |               |
%         |               |
% Neumann |               | Neumann
%         |               |
%         |               |
%       0 +---------------+------> x( i )
%         0   Dirichlet  2pi
%
% Spatial indexing: phi( i , j )
%   i: index along x: i = 1 , ... , Ni
%   j: index along y: j = 1 , ... , Ni
%
% Linear indexing: phi( k )
%   k = i + (j-1) * Ni
%
% The number of unknowns is Ni * Ni
% 
%   k_max = i_max + (j_max-1) * Ni =
%           Ni    + (Ni   -1) * Ni = Ni * Ni
%
% ------------------------------------------------------------------------

% Constant properties (rho, gamma) and advection velocity (w)
rho   = 1 ;
Gamma = 1 ;
w_x   = 0 ;
w_y   = 1 ;

% Domain square side length
L = 2*pi ;

% BC values
phi_BC = 0 ;
Jd_BC  = 0 ;

% FV number for each dimension and FV total number
Ni = 100 ;
N  = Ni * Ni ;

% FV side length
dx = L / Ni ;

% FV 'volume' (surface)
dA = dx * dx ;

% FV centroids coordinates
X = dx * ( (1:Ni) - 0.5 )' ;
Y = X ;

% FV equation:
%     A_P*phi_P + A_E*phi_E + A_W*phi_W + A_N*phi_N + A_S*phi_S = S_P
%
% Final system, in compact (matrix) form:
%                           A * phi = S
%
%   A: sparse coefficient matrix
% phi: solution vector
%   S: source term vector (RHS)

% FV equation coefficients (constants for each FV)
A_N =  rho * dx * w_y / 2 - Gamma ;
A_S = -rho * dx * w_y / 2 - Gamma ;
A_E =  rho * dx * w_x / 2 - Gamma ;
A_W = -rho * dx * w_x / 2 - Gamma ;
A_P = -( A_E + A_W + A_N + A_S ) ;

% Preparation of the 5 diagonals of A, now stored as Ni x Ni matrices
% (spatial indexing)
D_W = A_W * ones( Ni , Ni ) ;
D_E = A_E * ones( Ni , Ni ) ;
D_N = A_N * ones( Ni , Ni ) ;
D_S = A_S * ones( Ni , Ni ) ;
D_P = A_P * ones( Ni , Ni ) ;

% RHS preparation (midpoint second order integration)
% stored as Ni x Ni matrix (spatial indexing)
S = sin( X ) * sin( Y' ) * dA ;

% Diagonals correction in order to impose BC
%   South side ( y=0, j=1 ) (Dirichlet)
D_P( : , 1 ) = D_P( : , 1 ) - D_S( : , 1 ) ;
  S( : , 1 ) =   S( : , 1 ) - 2 * D_S( : , 1 ) * phi_BC ;
D_S( : , 1 ) = 0 ;

%   North side ( y=L, j=Ni ) (Neumann)
D_P( : , Ni ) = D_P( : , Ni ) + D_N( : , Ni ) ;
  S( : , Ni ) =   S( : , Ni ) - D_N( : , Ni ) * dx * Jd_BC / Gamma ;
D_N( : , Ni ) = 0 ;

%   West side ( x=0, i=1 ) (Neumann)
D_P( 1 , : ) = D_P( 1 , : ) + D_W( 1 , : ) ;
  S( 1 , : ) =   S( 1 , : ) - D_W( 1 , : ) * dx * Jd_BC / Gamma ;
D_W( 1 , : ) = 0 ;

%   East side ( x=L, i=Ni ) (Neumann)
D_P( Ni , : ) = D_P( Ni , : ) + D_E( Ni , : ) ;
  S( Ni , : ) =   S( Ni , : ) - D_E( Ni , : ) * dx * Jd_BC / Gamma ;
D_E( Ni , : ) = 0 ;

% Diagonals (and RHS) reshape to get vectors (linear indexing)
D_P = reshape( D_P , N , 1 ) ;
D_W = reshape( D_W , N , 1 ) ;
D_E = reshape( D_E , N , 1 ) ;
D_N = reshape( D_N , N , 1 ) ;
D_S = reshape( D_S , N , 1 ) ;
S   = reshape(   S , N , 1 ) ;

% Lower diagonals (South/West) and upper diagonals (East/North) shift.
% This shift operation is needed because of the way the following spdiags()
% command operates.
% (The circular property is not strictly needed here...)
D_S = circshift( D_S , -Ni ) ;
D_W = circshift( D_W ,  -1 ) ;
D_E = circshift( D_E ,   1 ) ;
D_N = circshift( D_N ,  Ni ) ;

% Coefficient matrix A creation from diagonals (sparse form)
A = spdiags( [ D_S D_W D_P D_E D_N ] , [ -Ni -1 0 1 Ni ] , N , N ) ;

% Direct solution (\ backslash)
phi = A\S ;

% Conjugate Gradient solution
% Coefficient matrix must be Symmetric Positive Definite (SPD).
% Be careful to the sign of A:
%     phi = pcg( A , S ) ;
%
% The use of a preconditioner heavily affects the rate of convergence; for
% example we can use an incomplete Cholesky factorization as
% preconditioner:
%     tol    = relative tolerance on residual
%     n_iter = maximum number of iterations
%     L = ichol( A ) ;
%     phi_pcg = pcg( A , S , tol , n_iter , L , L' ) ;

% Solution reshape to get again the spatial indexing
phi = reshape( phi , [ Ni , Ni ] ) ;

% Contour line plot (orthogonal to the diffusive flux) on figure 1
figure( 1 ) ;
contourf( X , Y , phi' ) ;

% Naming axes
xlabel( 'x' ) ;
ylabel( 'y' ) ;

% Surface plot on figure 2
figure( 2 ) ;
Superf = surf( X , Y , phi' ) ;

% Naming axes
xlabel( 'x' ) ;
ylabel( 'y' ) ;
zlabel( 'phi' ) ;

% Pretty surface properties
Superf.EdgeColor    = 'none'    ; % FV sides not drawn
Superf.FaceLighting = 'gouraud' ; % Lighting model
Superf.FaceColor    = 'interp'  ; % Color interpolation inside FV

% Adding a light
Light_1          = light ;
Light_1.Position = [ -2 0 2 ] ;
Light_1.Color    = [ .7 .7 .7 ] ;