Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/matlab/test_transf2d_cyl.m
c723281d03b6feature/refactoring

Last Change
Raw File

test_transf2d_cyl.m
View Options

%
% @file test_transf2d_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/>.
%
% @author
% (in alphabetical order)
% @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
%
file='test_transf2d_cyl.h5';
%
nx(1)=h5readatt(file,'/','NX');
ny(1)=h5readatt(file,'/','NY');
modem=h5readatt(file,'/','MODEM');
modep=h5readatt(file,'/','MODEP');
nidbas1=h5readatt(file,'/','NIDBAS1');
nidbas2=h5readatt(file,'/','NIDBAS2');
nlevels=h5readatt(file,'/','LEVELS');
title_str=sprintf('N=(%d,%d), NIDBAS=(%d,%d), MODEM=%d, MODEP=%d', ...
nx(1),ny(1),nidbas1,nidbas2,modem,modep);
%
% Grid sizes on each levels
%
for l=2:nlevels
nx(l) = nx(l-1)/2;
ny(l) = ny(l-1)/2;
end
%
% 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);
x=h5read(file,strcat(mglevel,'/x'));
y=h5read(file,strcat(mglevel,'/y'));
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
%
sol_direct=h5read(file,'/solutions/direct');
sol_anal=h5read(file,'/solutions/anal');
vfine=h5read(file,'/solutions/vfine');
figure
subplot(211)
surf(x,y,sol_direct')
xlabel('X'); ylabel('Y');
title('Direct Solve on the finest grid')
subplot(212)
surf(x,y,vfine')
xlabel('X'); ylabel('Y');
title('Prolongation solution on the finest grid')
figure
subplot(311)
plot(x, sol_direct(:,ny(1)/2),x, vfine(:,ny(1)/2),'o')
xlabel('r');
grid on
legend('Direct Solution', 'Prolonged Solution')
title(title_str)
subplot(313)
plot(y, sol_direct(1,:),y, vfine(1,:),'o')
xlabel('\theta'); ylabel('On axis')
grid on
title(title_str)
subplot(312)
plot(y, sol_direct(nx(1)/2,:),y, vfine(nx(1)/2,:),'o')
xlabel('\theta'); ylabel('Off axis')
grid on
title(title_str)
%
% Errors
%
errdisc = h5read(file,'/errors/errdisc');
resid = h5read(file,'/errors/resid');
restrict = h5read(file,'/errors/restrict');
prolong = h5read(file,'/errors/prolong');
errdisc_prolong = h5read(file,'/errors/disc_err_prolong');
figure
loglog(nx, errdisc, 'o-', nx(2:end), errdisc_prolong, 'h-')
grid on;
legend('Direct Solution', 'Prolonged Solution')
xlabel('N'); ylabel('Discretization Errors')
title(title_str);
c4science · Help