Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/matlab/test_mg2d_cyl.m
03540012fedefeature/fft-split

Last Change
Raw File

test_mg2d_cyl.m
View Options

%
% @file test_mg2d_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_mg2d_cyl.h5';
%
nx=h5readatt(file,'/','NX');
ny=h5readatt(file,'/','NY');
nidbas1=h5readatt(file,'/','NIDBAS1');
nidbas2=h5readatt(file,'/','NIDBAS2');
modem=h5readatt(file,'/','MODEM');
modep=h5readatt(file,'/','MODEP');
relax=h5readatt(file,'/','RELAX');
nlevels=h5readatt(file,'/','LEVELS');
nu1=h5readatt(file,'/','NU1');
nu2=h5readatt(file,'/','NU2');
mu=h5readatt(file,'/','MU');
title_str=sprintf(['N=(%d,%d), NIDBAS=(%d,%d), relax=%s, nu1=%d, ' ...
'nu2=%d, mu=%d, LEVELS=%d, MODEM=%d, MODEP=%d'], ...
nx,ny,nidbas1,nidbas2,relax,nu1,nu2,mu,nlevels,modem,modep);
%
% Prolongation matrices at the coarsest grid
%
levels=nlevels;
mglevel=sprintf('/mglevels/level.%.2d', levels);
dset=strcat(mglevel,'/matpx');
matpx=csr_mat(file,dset);
dset=strcat(mglevel,'/matpy');
matpy=csr_mat(file,dset);
%
% FE matrix at the finest grid
%
levels=1;
mglevel=sprintf('/mglevels/level.%.2d', levels);
dset=strcat(mglevel,'/mata');
[mata,diag]=csr_mat(file,dset);
f=h5read(file,strcat(mglevel,'/f'));
v=h5read(file,strcat(mglevel,'/v'));
f1d=h5read(file,strcat(mglevel,'/f1d'));
v1d=h5read(file,strcat(mglevel,'/v1d'));
% $$$ figure
% $$$ spy(mata)
%
% Solutions at the finest grid
%
x=h5read(file,'/solutions/xg');
y=h5read(file,'/solutions/yg');
sol_anal=h5read(file,'/solutions/anal');
sol_calc=h5read(file,'/solutions/calc');
sol_direct=h5read(file,'/solutions/direct');
nx=int32(size(x,1));
ny=int32(size(y,1));
figure
subplot(211)
surf(x,y,sol_anal')
xlabel('X'); ylabel('Y');
title('Analytical solution on the finest grid')
subplot(212)
surf(x,y,sol_calc')
% $$$ surf(x,y, (abs(sol_calc'-sol_anal')))
xlabel('X'); ylabel('Y');
title('Calculated solution on the finest grid')
figure
subplot(211)
plot(x, sol_anal(:,ny/2),x, sol_calc(:,ny/2),'o')
xlabel('r');
grid on
legend('Analytic Solution', 'MG Solution')
title(title_str)
subplot(212)
if modem == 0
plot(y, sol_anal(1,:),y, sol_calc(1,:),'o')
else
plot(y, sol_anal(nx/2,:),y, sol_calc(nx/2,:),'o')
end
xlabel('\theta');
grid on
title(title_str)
% $$$ figure
% $$$ subplot(211)
% $$$ semilogy(x, abs(sol_anal(:,ny/2)-sol_calc(:,ny/2)),'o')
% $$$ xlabel('r'); ylabel('Error')
% $$$ grid on
% $$$ title(title_str)
% $$$ subplot(212)
% $$$ if modem == 0
% $$$ semilogy(y, abs(sol_anal(1,:)-sol_calc(1,:)),'o')
% $$$ else
% $$$ semilogy(y, abs(sol_anal(nx/2,:)-sol_calc(nx/2,:)),'o')
% $$$ end
% $$$ xlabel('\theta'); ylabel('Error');
% $$$ grid on
% $$$ title(title_str)
% $$$
%
% Iterations
%
dset='/Iterations/';
disc_err=h5read(file, strcat(dset,'disc_errors'));
resid=h5read(file, strcat(dset,'residues'));
its=0:1:size(resid,1)-1;
figure
subplot(211)
semilogy(its,resid,'o-')
grid on
xlabel('Iterations'); ylabel('Norm of residue');
title(title_str);
subplot(212)
semilogy(its,disc_err,'h-')
grid on
xlabel('Iterations'); ylabel('Norm of error');
c4science · Help