%
% @file ppoisson_fd.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>
%
if ~exist('file'), file='ppoisson_fd.h5'; end

prb=h5readatt(file,'/','PRB');
nx=h5readatt(file,'/','NX');
ny=h5readatt(file,'/','NY');
lx=h5readatt(file,'/','LX');
ly=h5readatt(file,'/','LY');
kx=h5readatt(file,'/','KX');
ky=h5readatt(file,'/','KY');
beta=h5readatt(file,'/','BETA');
omega=h5readatt(file,'/','OMEGA');
relax=h5readatt(file,'/','RELAX');
nlevels=h5readatt(file,'/','LEVELS');
mu=h5readatt(file,'/','MU');
nu1=h5readatt(file,'/','NU1');
nu2=h5readatt(file,'/','NU2');
direct_solve_nits=h5readatt(file,'/','DIRECT_SOLVE_NITS');

title_str=sprintf(['PRB=%s, N=(%d,%d), relax=%s, V(%d,%d), LEVELS=%d, DIRECT SOLVE=%d'], ...
                  prb, nx, ny, relax, nu1, nu2, nlevels, direct_solve_nits);

x = h5read(file, '/xgrid');
y = h5read(file, '/ygrid');
[X,Y]=meshgrid(y,x);
n1=size(x,1);
n2=size(y,1);
n=n1*n2

mat = stencil_mat(file, '/MAT');
% $$$ figure
% $$$ spy(mat)

f = h5read(file,'/f'); f1d=reshape(f,n,1);
v = h5read(file,'/v'); v1d=reshape(v,n,1);
u = h5read(file,'/u'); u1d=reshape(u,n,1);

% $$$ udirect1d = mat\f1d;
% $$$ udirect=reshape(udirect1d,n1,n2);
% $$$ fprintf('Residual of direct solution = %.3e\n', norm(mat*udirect1d-f1d));
% $$$ fprintf('Error of direct solution    = %.3e\n', norm(udirect1d- ...
% $$$                                                   v1d));
figure
subplot(221)
pcolor(x,y,v'-u')
xlabel('X'); ylabel('Y');
shading interp
colorbar
title('Error on the finest grid')

subplot(222)
[yy,iy] = max(abs(v),[],2);
[xx,ix] = max(yy);
iy0=iy(ix);
str=sprintf('Solution at y = %.4f', y(iy0));
plot(x, v(:,iy0),x, u(:,iy0),'o')
xlabel('x'); ylabel(str);
grid on
legend('Analytic Solution', 'Computed Solution')
title(title_str)

subplot(223)
[xx,ix] = max(abs(v));
[yy,iy] = max(xx);
ix0=ix(iy);
str=sprintf('Solution at x = %.4f', x(ix0));
plot(y, v(ix0,:),y, u(ix0,:),'o')
xlabel('y'); ylabel(str);
grid on

figure
resid_it=h5read(file, '/resid');
err_it=h5read(file, '/error');
nits = size(resid_it,1)-1;
it=0:nits;
subplot(211)
semilogy(it,resid_it,'o-')
xlabel('Iterations')
ylabel('Residual norm')
grid on
title(title_str)
subplot(212)
semilogy(it, err_it, 'o-')
xlabel('Iterations')
ylabel('Norm of Discretization Error')
grid on