
%
% Apply preconditioned GMRES to solve the discretized advection diffusion equation
%
%     - diff*Delta y(x) + adv'* nabla y(x) + react y(x)  = f(x)  x in (0,1)^2,   
%                                                  y(x)  = 0     x on the boundary
%          
% We use f = 1.
%
% Matthias Heinkenschloss
% March 13, 2009 (last modified February 7, 2017)



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



% parameters in the advection diffusion equation
diff   = 1.e-2;   % must be positive
adv    = [1; 1];  % 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)';

%% 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 + 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 + 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);



% right hand side
b = ones(n1*n2,1);

    
% FD  solution
yfd = A\b;

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

figure(1)
surf( x1, x2, Y'); shading interp;
xlabel('x_1'); ylabel('x_2'); view(-60,45);

    


% Get matrics for pointwise preconditioners
E = tril(-A,-1);
F = triu(-A,1);
D = spdiags(diag(A),0,n1*n2,n1*n2);

%% GMRES options
options.tol     = 1.e-7;
options.maxit   = 200;
options.restart = 201;

y0 = zeros(size(b));

%% solve the system using GMRES with left precondtioning
% solve using GMRES without preconditioning
[y, res, iter, flag] = mygmres( A, b, [], [], y0, options );

figure(2)
fprintf(' GMRES (no prec.) returned after %d iterations with flag = %d\n', iter,flag)
semilogy((0:iter),res,'k-'); hold on

% solve using GMRES with left Jacobi preconditioning
[y, res, iter, flag] = mygmres( A, b, D, [], y0, options );

figure(2)
fprintf(' GMRES (left Jacobi prec.) returned after %d iterations with flag = %d\n', iter,flag)
semilogy((0:iter),res,'r--'); hold on


% solve using GMRES with left GS preconditioning
M     = D-E;

[y, res, iter, flag] = mygmres( A, b, M, [], y0, options );

figure(2)
fprintf(' GMRES (left GS prec.) returned after %d iterations with flag = %d\n', iter,flag)
semilogy((0:iter),res,'b-.'); hold on

% solve using GMRES with left SOR preconditioning
omega = 1.2;
M     = ((1/omega)*D-E);

[y, res, iter, flag] = mygmres( A, b, M, [], y0, options );

figure(2)
fprintf(' GMRES (left SOR prec.) returned after %d iterations with flag = %d\n', iter,flag)
semilogy((0:iter),res,'g-'); hold off

legend('no prec.', 'left Jacobi prec.', 'left GS prec.', 'left SOR prec.',...
       'Location','SouthEast')
axis('tight')
%print -depsc advdiff2D_left_prec_gmres_res.eps


%% solve the system using GMRES with right precondtioning
figure(3)
% solve using GMRES without preconditioning
[x, res, iter, flag] = mygmres( A, b, [], [], y0, options );



fprintf(' GMRES (no prec.) returned after %d iterations with flag = %d\n', iter,flag)
semilogy((0:iter),res,'k-'); hold on

% solve using GMRES with right Jacobi preconditioning
[y, res, iter, flag] = mygmres( A, b, [], D, y0, options );


fprintf(' GMRES (right Jacobi prec.) returned after %d iterations with flag = %d\n', iter,flag)
semilogy((0:iter),res,'r--'); hold on


% solve using GMRES with right GS preconditioning
M     = D-E;

[y, res, iter, flag] = mygmres( A, b, [], M, y0, options );


fprintf(' GMRES (right GS prec.) returned after %d iterations with flag = %d\n', iter,flag)
semilogy((0:iter),res,'b-.'); hold on

% solve using GMRES with right SOR preconditioning
omega = 1.2;
M     = ((1/omega)*D-E);

[y, res, iter, flag] = mygmres( A, b, [], M, y0, options );


fprintf(' GMRES (right SOR prec.) returned after %d iterations with flag = %d\n', iter,flag)
semilogy((0:iter),res,'g-'); hold off

legend('no prec.', 'right Jacobi prec.', 'right GS prec.', 'right SOR prec.',...
       'Location','SouthEast')
axis('tight')

%print -depsc advdiff2D_right_prec_gmres_res.eps

