% ------------------------------------------------------------------------
% 
% Problem: 1D steady-state advection-diffusion, harmonic RHS:
%                (rho*w*phi)_x = (Gamma*phi_x)_x + sin(x)
%
% Domain: [0,2pi]
%
% Grid: uniform (costant spacing)
%
% Boundary conditions (BC):
%   Dirichlet ( phi           = phi_BC ) on left  end ( x = 0   )
%   Neumann   ( Gamma * phi_x = Jd_BC  ) on right end ( x = 2pi )
%
% Discretization method: Finite Volume (FV)
%
% Solution methodology: direct (sparse coefficient matrix)
%
% Date: 30/3/2016
% Author: R. Zamolo
%
% ------------------------------------------------------------------------
%
% Domain, coordinates systems and boundary conditions:
%
% Dirichlet        Neumann
%   --|---------------|------> x( i )
%     0              2pi
%
% Spatial indexing: phi( i )
%   i: index along x: i = 1 , ... , N
%
% ------------------------------------------------------------------------

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

% Domain length
L = 2*pi ;

% BC values
phi_BC = 0 ;
Jd_BC  = 0 ;

% FV number
N = 100 ;

% FV side length
dx = L / N ;

% FV centroids abscissae (column vector)
X = dx * ( (1:N) - 0.5 )' ;

% FV equation:
%              A_W*phi_W + A_P*phi_P + A_E*phi_E = 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_W = -rho * w / 2 - Gamma / dx ;
A_E =  rho * w / 2 - Gamma / dx ;
A_P = -( A_W + A_E ) ;

% Preparation of the 3 diagonals of A, now stored as (column) vectors
D_W = A_W * ones( N , 1 ) ;
D_P = A_P * ones( N , 1 ) ;
D_E = A_E * ones( N , 1 ) ;

% RHS vector preparation (midpoint second order integration)
S = sin( X ) * dx ;

% Diagonals correction in order to impose BC
%   Left  end ( x=0, i=1 ) (Dirichlet)
D_P( 1 ) = D_P( 1 ) - D_W( 1 ) ;
  S( 1 ) =   S( 1 ) - 2 * D_W( 1 ) * phi_BC ;
D_W( 1 ) = 0 ;

%   Right end ( x=L, i=N ) (Neumann)
D_P( N ) = D_P( N ) + D_E( N ) ;
  S( N ) =   S( N ) - D_E( N ) * dx * Jd_BC / Gamma ;
D_E( N ) = 0 ;

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

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

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

% Solution plot
plot( X , phi ) ;

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