Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/matlab/test_jacobi.m
c723281d03b6feature/refactoring

Last Change
Raw File

test_jacobi.m
View Options

%
% @file test_jacobi.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_jacobi.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);
resids = h5read(file,'/resids');
errs = h5read(file,'/errs');
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(211)
[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(212)
[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
% $$$ subplot(321)
% $$$ surf(x,y,v'); title('Exact solution')
% $$$ xlabel('X'); ylabel('Y')
% $$$ subplot(322)
% $$$ surf(x,y,f'); title('RHS')
% $$$ xlabel('X'); ylabel('Y')
% $$$ subplot(323)
% $$$ surf(x,y,u'); title('Num. solution')
% $$$ subplot(324)
% $$$ surf(x,y,udirect'); title('Direct. solution')
% $$$ xlabel('X'); ylabel('Y')
% $$$ subplot(325)
% $$$ surf(x,y,resids'); title('Residuals')
% $$$ xlabel('X'); ylabel('Y')
% $$$ subplot(326)
% $$$ surf(x,y,errs'); title('Errors')
% $$$ xlabel('X'); ylabel('Y')
resid_it=h5read(file, '/resid');
err_it=h5read(file, '/error');
nits = size(resid_it,1)-1;
it=0:nits;
figure
semilogy(it,resid_it, it, err_it)
legend('Residual norm', 'Discretization error')
xlabel('Iterations')
grid on
c4science · Help