% ------------------------------------------------------------------------
% 
% Problem: 2D Poisson, harmonic Right Hand Side (RHS):
%                   phi_xx + phi_yy = cos(x) * cos(y)
%
% Domain: square [0,2pi] x [0,2pi]
%
% Grid: structured, cartesian, constant spacing, dx = dy
%
% Boundary conditions (BC):
%   Neumann ( phi_n = Dphi/Dn = 0 )
%   phi(0,0) = 0 (required to get unique solution)
%
% Discretization method: Finite Volume (FV)
%
% Solution methodology: iterative Jacobi/
%                                 SOR (Successive Over-Relaxation)/
%                                 Red-Black SOR
%
% Date: 31/3/2016
% Author: R. Zamolo
%
% ------------------------------------------------------------------------
%
% Domain, coordinates systems and boundary conditions:
%
%         ^ y( j ) 
%         |
%         |
%         |    phi_n=0
%     2pi +---------------+
%         |               |
%         |               |
% phi_n=0 |               | phi_n=0
%         |  phi=0        |
%         | /             |
%       0 +---------------+------> x( i )
%         0    phi_n=0   2pi
%
% Spatial indexing: phi( i , j )
%   i: index along x: i = 1 , ... , Ni
%   j: index along y: j = 1 , ... , Ni
%
% ------------------------------------------------------------------------

% Domain square side length
L = 2*pi ;

% FV number for each dimension
Ni = 50 ;

% 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

% FV equation coefficients (constants for each FV)
A_W =  1 ;
A_E =  1 ;
A_N =  1 ;
A_S =  1 ;
A_P = -4 ;

% Solution initialization (spatial indexing)
% The values of the 'ghost cells' unknowns are required to be stored, so
% the solution matrix phi is (Ni+2) x (Ni+2)
phi = zeros( Ni+2 , Ni+2 ) ;

% Source term Ni x Ni matrix (spatial indexing)
S = cos( X ) * cos( Y' ) * dA ;

% Indexes of FV inside the domain (not ghost cells) => phi(k,k) 
k = 2 : (Ni+1) ;

% Number of iterations
n_iter = 1000 ;

% SOR/Red_Black over-relaxation parameter
w0 = 1.5 ; % <2

% Jacobi=1 / SOR=2 / Red-Black SOR=3, selector
Method = 3 ;

% Indexes for Red-Black Ordering
si = [ 0 1 0 1 ] ;
sj = [ 0 1 1 0 ] ;

% Residual interval
d_iter_Res = 8 ;

% Residual RMS
err = zeros( n_iter / d_iter_Res , 1 ) ;

% Iterative cycle
tic() ;
for iter = 1 : n_iter

  % Multi w
  w = w0 ;
  %if mod( iter , 80 )==0 && Method==3 ; w = 5 ; end % RB-SOR
  
  % Updating FV values
  switch Method
    case 1
      % Jacobi iteration
      ngbrs = A_W*phi( k-1 , k ) + A_E*phi( k+1 , k ) + ...
              A_N*phi( k , k+1 ) + A_S*phi( k , k-1 ) ;
      phi( k , k ) = ( S - ngbrs ) / A_P ;
    case 2
      % SOR iteration
      for i = k
        for j = k
          ngbrs = A_W*phi( i-1 , j ) + A_E*phi( i+1 , j ) + ...
                  A_N*phi( i , j+1 ) + A_S*phi( i , j-1 ) ;
          phi( i , j ) = w*( S(i-1,j-1) - ngbrs )/A_P + (1-w)*phi( i , j ) ;
        end
      end
    case 3
      % Red-Black SOR iteration (Red cells: RB=1,2; Black cells: RB=3,4)
      for RB = 1 : 4
        i = (2+si(RB)) : 2 : (Ni+1) ;
        j = (2+sj(RB)) : 2 : (Ni+1) ;
        ngbrs = A_W*phi( i-1 , j ) + A_E*phi( i+1 , j ) + ...
                A_N*phi( i , j+1 ) + A_S*phi( i , j-1 ) ;
        phi( i , j ) = w*( S(i-1,j-1) - ngbrs ) / A_P + (1-w)*phi( i , j ) ;
      end
  end
  
  % Updating 'ghost cells' values to enforce BC
  %   South side ( y=0, j=1 ) (Neumann)
  phi( k , 1 ) = phi( k , 2 ) ;

  %   North side ( y=L, j=end ) (Neumann)
  phi( k , end ) = phi( k , end-1 ) ;
  
  %   West side ( x=0, i=1 ) (Neumann)
  phi( 1 , k ) = phi( 2 , k ) ;

  %   East side ( x=L, i=end ) (Neumann)
  phi( end , k ) = phi( end-1 , k ) ;
  
  % Null solution on first FV
  phi = phi - phi(2,2) ;
  
  % Residual
  if mod( iter , d_iter_Res ) == 0
    Res = A_W*phi( k-1 , k ) + A_E*phi( k+1 , k ) + ...
          A_N*phi( k , k+1 ) + A_S*phi( k , k-1 ) + A_P*phi( k , k ) - S ;
    err( iter / d_iter_Res ) = norm( Res , 'fro' ) / norm( S , 'fro' ) ;
  end
  
end
toc() ;

% Residual history plot
% on figure 1
figure( 1 ) ;
hold on ;
plot( 1:d_iter_Res:n_iter , err ) ;
ax = gca() ;
ax.YScale = 'log' ;
hold on ;
fig = gcf ;
fig.Position = [ 10 350 560 420 ] ;

% Naming axes
xlabel( 'Iterations' ) ;
ylabel( 'Residual error (RMS)' ) ;


% Extraction of solution values for FV inside the domain
phi = phi( k , k ) ;

% Contour line plot (orthogonal to the pressure gradient)
% on figure 2
figure( 2 ) ;
contourf( X , Y , phi' ) ;

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

fig = gcf ;
fig.Position = [ 900 300 560 420 ] ;

% Surface plot
% on figure 3
figure( 3 ) ;
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 ] ;

fig = gcf ;
fig.Position = [ 900 50 560 420 ] ;