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

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')

subplot(223)
[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')

subplot(224)
[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

resid_it=h5read(file, '/resid');
err_it=h5read(file, '/error');
nits = size(resid_it,1)-1;
it=0:nits;

subplot(222)
semilogy(it,resid_it, it, err_it)
legend('Residual norm', 'Discretization error')
xlabel('Iterations')
grid on