%
% @file test_mgp.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/>.
%
% @author
% (in alphabetical order)
% @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
%
file='test_mgp.h5';
%
nx=h5readatt(file,'/','NX');
nidbas=h5readatt(file,'/','NIDBAS');
relax=h5readatt(file,'/','RELAX');
levels=h5readatt(file,'/','LEVELS');
nu1=h5readatt(file,'/','NU1');
nu2=h5readatt(file,'/','NU2');
title_str=sprintf('NX = %d, NIDBAS = %d, levels = %d, nu1 = %d, nu2 = %d', nx, nidbas, levels, nu1, nu2);
%
% Read matrices at coarset grid
%
for lev=2:levels
%
%  FE mat at fine grid
mglevel=sprintf('/mglevels/level.%.2d', lev-1);
dset=strcat(mglevel,'/mata');
mata_f=h5read(file,dset);
n=size(mata_f,1);
%
%  FE mat at coarse grid
mglevel=sprintf('/mglevels/level.%.2d', lev);
dset=strcat(mglevel,'/mata');
mata_c=h5read(file,dset);
n=size(mata_c,1);
%
%  Prolong mat
dset=strcat(mglevel,'/matp');
matp=h5read(file,dset);
%
%   Check
fprintf(1,'Level %d: ||A_coarse - P''*A_fine*P|| = %g\n', lev, norm(matp'*mata_f*matp ...
                                                  - mata_c))
end
%
%   Iterations
dset='/Iterations/';
err=h5read(file, strcat(dset,'errors'));
disc_err=h5read(file, strcat(dset,'disc_errors'));
resid=h5read(file, strcat(dset,'residues'));
its=0:1:size(err,1)-1;

figure
subplot(221)
semilogy(its,resid,'o-', its, disc_err,'h-')
legend('Residue', 'Error')
grid on
xlabel('Iterations'); ylabel('Norm od residue and error');
title(title_str);

%
%   Plot grid values at the last iteration
xgrid=h5read(file, '/Iterations/xgrid');
u_calc=h5read(file, '/Iterations/u_calc');
u_exact=h5read(file, '/Iterations/u_exact');
u_direct=h5read(file, '/Iterations/u_direct');
subplot(222)
plot(xgrid,u_exact, xgrid,u_calc,'o')
xlabel('X');ylabel('Grid values of solution')
grid on
title(title_str);

subplot(223)
semilogy(xgrid,abs(u_calc-u_direct))
xlabel('X');ylabel('Diff with direct solution')
grid on
title(title_str);

subplot(224)
semilogy(xgrid,abs(u_calc-u_exact))
xlabel('X');ylabel('Diff with exact solution')
grid on
title(title_str);