%
% Finite difference approximation to the solution of the 2D advection diffusion equation
%
%     - diff Delta y(x) + adv . nabla y(x) + react y(x) = f(x)  x in (0,1)^2,   
%                                                  y(x) = g(x)  x on the boundary
%          
%
% Matthias Heinkenschloss
% Jan 5, 2017
%

set(0, 'defaultaxesfontsize',24,'defaultaxeslinewidth',2.0,...
       'defaultlinelinewidth',2.0,'defaultpatchlinewidth',2.0,...
       'defaulttextfontsize',24);


clear all

diff   = 1.e-4;                    % must be positive
theta  = 47.3 / 180 * pi;
adv    = [cos(theta),sin(theta)];  % must be component wise nonnegative
%    (the upwind scheme is coded for component wise nonnegative advection)
react  = 0;                        % must be nonnegative

% set mesh points
n1  = 100;
n2  = 100;
h1  = 1/(n1+1);
h2  = 1/(n2+1);
x1  = (0:h1:1)';
x2  = (0:h2:1)';

kronecker = 2;  % = 1 use Kronecker structure; = 0 use blockdiagonal
if(kronecker)
    %% Generate system matrix (lexicographic ordering)
    %  using Kronecker products
    %  Generate A1
    c    = -(diff/h1^2 + adv(1)/h1)*ones(n1,1);
    d    = (2*diff/h1^2 + adv(1)/h1 + 0.5*react)*ones(n1,1);
    e    = (-diff/h1^2)*ones(n1,1);
    A1   = spdiags([c d e], -1:1, n1, n1);
    
    %  Generate A2
    c    = -(diff/h2^2 + adv(2)/h2)*ones(n2,1);
    d    = (2*diff/h2^2 + adv(2)/h2 + 0.5*react)*ones(n2,1);
    e    = (-diff/h2^2)*ones(n2,1);
    A2   = spdiags([c d e], -1:1, n2, n2);
    
    I1   = speye(n1,n1);
    I2   = speye(n2,n2);
    A    = kron(I2,A1) + kron(A2,I1);
    
else
    %% Generate system matrix (lexicographic ordering)
    %  explicityly as a block tridiagonal matrix
    % set subdiagonal blocks
    e     = ((diff+h2*adv(2))/h2^2)*ones(n1,1);
    Eb    = spdiags(e,0,n1,n1);
    % set diagonal blocks
    e1    = (-(diff+h1*adv(1))/h1^2)*ones(n1,1);
    e2    = ((2*diff+h1*adv(1))/h1^2 + (2*diff+h2*adv(2))/h2^2 + react)*ones(n1,1);
    e3    = (-diff/h1^2)*ones(n1,1);
    Db    = spdiags([e1 e2 e3], -1:1, n1, n1);
    % set superdiagonal blocks
    e     = (diff/h2^2)*ones(n1,1);
    Fb    = spdiags(e,0,n1,n1);
    
    % generate system matrix
    A     = blktridiag(Db,-Eb,-Fb,n2);
end

figure(1)
spy(A)
title('Sparsity Pattern of A')


%% right hand side
% boundary data  g(i,j) = Dirichlet boundary data at (x_i, y_j)
g        = zeros(n1+2,n2+2);
g(:,1)   = 1;               % g(x,y) = 1 for y = 0
ind      = ( abs(x2) <= 0.25+eps );
g(1,ind) = 1;               % g(x,y) = 1 for x = 0 and y <= 0.25

% right hand side function  f(i,j) = right hand side function at (x_i, y_j)
f = zeros(n1+2,n2+2);

% right hand side of linear system
b                    = reshape( f(2:n1+1,2:n2+1), n1*n2, 1); 
b(1:n1)              = b(1:n1)              + (diff/h2^2 + adv(2)/h2)*g(2:n1+1,1);
b((n2-1)*n1+1:n2*n1) = b((n2-1)*n1+1:n2*n1) + (diff/h2^2)*g(2:n1+1,n2+2);
b(1:n1:n1*n2)        = b(1:n1:n1*n2)        + (diff/h1^2 + adv(1)/h1)*g(1,2:n2+1)';
b(n1:n1:n1*n2)       = b(n1:n1:n1*n2)       + (diff/h1^2)*g(n1+2,2:n2+1)';


%% Compute and plot the solution
yfd = A\b;   % 




% Y(i,j)  is the approximate solution at (x_i, y_j)
Y                = g;                      % set boundary data
Y(2:n1+1,2:n2+1) = reshape (yfd, n1, n2);  % solution in interior

figure(2)
surf( x1, x2, Y'); shading interp;
xlabel('x_1'); ylabel('x_2'); view(-60,45);
%print -depsc adv-diff-2d-ex1.eps
    

