Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/matlab/test_intergrid.m
93a1e9712d04feature/refactoring

Last Change
Raw File

test_intergrid.m
View Options

%
% @file test_intergrid.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='test_intergrid0.h5'; end
%
nx=h5readatt(file,'/','NX');
ny=h5readatt(file,'/','NY');
lx=h5readatt(file,'/','LX');
ly=h5readatt(file,'/','LY');
kx=h5readatt(file,'/','KX');
ky=h5readatt(file,'/','KY');
nlevels=h5readatt(file,'/','LEVELS');
title_str=sprintf(['N=(%d,%d), Lx=%d, Ly=%d, LEVELS=%d, KX=%d, KY=%d'], ...
nx,ny,lx,ly,nlevels,kx,ky);
if nlevels ~= 2
disp 'levels should be 2!'
return
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);
for l=1:2
mglevel=sprintf('/mglevels/level.%.2d', l);
x=h5read(file,strcat(mglevel,'/x'));
y=h5read(file,strcat(mglevel,'/y'));
f=h5read(file,strcat(mglevel,'/f'));
v=h5read(file,strcat(mglevel,'/v'));
figure
subplot(221)
[yy,iy] = max(abs(f),[],2); [xx,ix] = max(yy); iy0=iy(ix);
str=sprintf('f at y = %.4f', y(iy0));
plot(x, f(:,iy0),'o-')
xlabel('x'); ylabel(str);
grid on
title(title_str)
subplot(222)
[xx,ix] = max(abs(f)); [yy,iy] = max(xx); ix0=ix(iy);
str=sprintf('f at x = %.4f', x(ix0));
plot(y, f(ix0,:),'o-')
xlabel('y'); ylabel(str);
grid on
title(title_str)
subplot(223)
[yy,iy] = max(abs(v),[],2); [xx,ix] = max(yy); iy0=iy(ix);
str=sprintf('v at y = %.4f', y(iy0));
plot(x, v(:,iy0),'ro-')
xlabel('x'); ylabel(str);
grid on
title(title_str)
subplot(224)
[xx,ix] = max(abs(v)); [yy,iy] = max(xx); ix0=ix(iy);
str=sprintf('v at x = %.4f', x(ix0));
plot(y, v(ix0,:),'ro-')
xlabel('y'); ylabel(str);
grid on
title(title_str)
if l==1
ffine=f; vfine=v;
else
fcoarse=f; vcoarse=v;
end
end
%% Check
err_restriction = matpx'*ffine*matpy./4 - fcoarse;
err_prolong = matpx*vcoarse*matpy' - vfine;
max(max(abs(err_restriction)))
max(max(abs(err_prolong)))
c4science · Help