%
% @file test_relax2d_cyl.m
%
% @brief
%
% @copyright
% Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
% SPC (Swiss Plasma Center)
%
% spclibs is free software: you can redistribute it and/or modify it under
% the terms of the GNU Lesser General Public License as published by the Free
% Software Foundation, either version 3 of the License, or (at your option)
% any later version.
%
% spclibs is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this program. If not, see <https://www.gnu.org/licenses/>.
%
% @authors
% (in alphabetical order)
% @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
%
file='test_relax2d_cyl.h5';
%
nx=h5readatt(file,'/','NX');
ny=h5readatt(file,'/','NY');
modem=h5readatt(file,'/','MODEM');
modep=h5readatt(file,'/','MODEP');
nidbas1=h5readatt(file,'/','NIDBAS1');
nidbas2=h5readatt(file,'/','NIDBAS2');
levels=h5readatt(file,'/','LEVELS');
omega=h5readatt(file,'/','OMEGA');
relax=h5readatt(file,'/','RELAX');
if relax(1:2) == 'ja'
    title_str=sprintf('N=(%d,%d), NIDBAS=(%d,%d), MODEM=%d, MODEP=%d, relax=%s, omega=%.3f', ...
                  nx,ny,nidbas1,nidbas2,modem,modep,relax, omega);
else
    title_str=sprintf('N=(%d,%d), NIDBAS=(%d,%d), MODEM=%d, MODEP=%d, relax=%s,', ...
                  nx,ny,nidbas1,nidbas2,modem,modep,relax);
end

%
%   Prolongation matrices at the coarsest grid
%
mglevel=sprintf('/mglevels/level.%.2d', levels);
dset=strcat(mglevel,'/matpx');
matpx=csr_mat(file,dset);
dset=strcat(mglevel,'/matpy');
matpy=csr_mat(file,dset);

%
%    Solutions at the finest grid
%
mglevel=sprintf('/mglevels/level.%.2d', 1);
x=h5read(file,strcat(mglevel,'/x'));
y=h5read(file,strcat(mglevel,'/y'));
f=h5read(file,strcat(mglevel,'/f'));
v=h5read(file,strcat(mglevel,'/v'));
f1d=h5read(file,strcat(mglevel,'/f1d'));
v1d=h5read(file,strcat(mglevel,'/v1d'));
sol_direct=h5read(file,'/solutions/direct');
sol_relax=h5read(file,'/solutions/relax');
sol_anal=h5read(file,'/solutions/anal');

% $$$ figure
% $$$ surf(x,y,sol_direct')
% $$$ xlabel('r'); ylabel('\theta');
% $$$ title(title_str)

% $$$ figure
% $$$ subplot(211)
% $$$ plot(x, sol_anal(:,ny/2),x, sol_direct(:,ny/2),'o')
% $$$ xlabel('r'); ylabel('Direct solution')
% $$$ grid on
% $$$ title(title_str)
% $$$ subplot(212)
% $$$ if modem == 0
% $$$     plot(x, sol_anal(1,:),x, sol_direct(1,:),'o')
% $$$ else
% $$$     plot(x, sol_anal(nx/2,:),x, sol_direct(nx/2,:),'o')
% $$$ end
% $$$ xlabel('\theta'); ylabel('Direct solution')
% $$$ grid on
% $$$ title(title_str)

%
%    Relaxations
%
errdisc=h5read(file,'/relaxation/errdisc');
resid=h5read(file,'/relaxation/resid');
its=0:1:size(errdisc)-1;
figure
semilogy(its,errdisc,its,resid)
legend('Discretisation Error', 'Residual Norm')
xlabel('Iterations')
grid on
title(title_str)

% $$$ figure
% $$$ subplot(211)
% $$$ plot(x, sol_anal(:,ny/2),x, sol_relax(:,ny/2),'o')
% $$$ xlabel('r'); ; ylabel('Relaxed solution')
% $$$ grid on
% $$$ title(title_str)
% $$$ subplot(212)
% $$$ if modem == 0
% $$$     plot(x, sol_anal(1,:),x, sol_relax(1,:),'o')
% $$$ else
% $$$     plot(x, sol_anal(nx/2,:),x, sol_relax(nx/2,:),'o')
% $$$ end
% $$$ xlabel('\theta'); ylabel('Relaxed solution')
% $$$ grid on
% $$$ title(title_str)

figure
subplot(211)
plot(x, sol_anal(:,ny/2),x, sol_direct(:,ny/2),'o')
xlabel('r'); ; ylabel('Direct solution')
grid on
title(title_str)
subplot(212)
plot(x, sol_anal(1,:),x, sol_direct(1,:),'o')
xlabel('\theta'); ylabel('Direct solution at the axis')
grid on
title(title_str)

%
%   FE matrix at the finest grid
%
levels=1;
mglevel=sprintf('/mglevels/level.%.2d', levels);
dset=strcat(mglevel,'/mata');
[mata,diag]=csr_mat(file,dset);
n=size(diag,1);
% $$$ figure
% $$$ spy(mata)
% $$$ title(title_str)

if relax(1:2) == 'ja'
    %
    %   Eigenvalues of inv(D)*A
    %
    matd=spdiags(diag,0,n,n);
    lambda = eigs(mata, matd, n);
    om=0:0.01:0.8;
    clear rho;
    for i=1:size(om,2)
        rho(i) = max(abs(1-om(i).*lambda));
    end
    figure
    subplot(211)
    plot(lambda,'o')
    ylabel('Eigenvalue of inv(D)*A')
    grid on
    title(title_str)
    subplot(212)
    plot(om,rho)
    xlabel('\omega'); ylabel('\rho')
    grid on
    omega_c = 2.0/max(lambda);
    lambda_min = eigs(mata, matd, 1, 'SM');
    fprintf(1, 'Spectral Radius of relaxation matrix = %g\n', abs(1-omega*lambda_min))
    str = sprintf('Critical omega = %.3f', omega_c)
    title(str)
else
%
%    Spectral radius of GS Iteration Matrix
%         Rg = -(D+L)^(-1) * U
%
    matl= tril(mata,0); % D+L
    lambda = eigs(-triu(mata,1),matl);
    fprintf(1, 'Spectral Radius of relaxation matrix = %g\n', max(abs(lambda)))
    figure
    plot(lambda, 'o')  
    axis equal
    grid on
    xlabel('Real of eigenvalues'); ylabel('Imag of eigenvalues')
    title(title_str)
end