Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/matlab/poisson_mg.m
93a1e9712d04feature/refactoring

Last Change
Raw File

poisson_mg.m
View Options

%
% @file poisson_mg.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='poisson_mg.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');
beta=h5readatt(file,'/','BETA');
omega=h5readatt(file,'/','OMEGA');
relax=h5readatt(file,'/','RELAX');
mat_type=h5readatt(file,'/','MAT_TYPE');
nlevels=h5readatt(file,'/','LEVELS');
mu=h5readatt(file,'/','MU');
nnu=h5readatt(file,'/','NNU');
nu1=h5read(file,'/nu1');
nu2=h5read(file,'/nu2');
title_str=sprintf(['N=(%d,%d), Lx=%d, Ly=%d, beta=%.4f, relax=%s, V(%d,%d), ' ...
'LEVELS=%d, KX=%d, KY=%d'], ...
nx,ny,lx,ly,beta,relax,nu1(nnu),nu2(nnu),nlevels,kx,ky);
%
% 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');
if mat_type == 'csr'
[mata,diag]=csr_mat(file,dset);
else
[mata,diag]=cds_mat(file,dset);
end
n=size(diag,1);
% $$$ figure
% $$$ spy(mata)
%
% Spectral radius of GS Iteration Matrix
% Rg = -(D+L)^(-1) * U
%
% $$$ matl= tril(mata,0); % D+L
% $$$ lambda = eigs(-triu(mata,1),matl)
% $$$ fprintf(1, 'Spectral Radius of GS relaxation matrix = %g\n', max(abs(lambda)))
% $$$ figure
% $$$ plot(lambda, 'o')
% $$$ axis equal
% $$$ grid on
% $$$ xlabel('Real of eigenvalues'); ylabel('Imag of eigenvalues')
% $$$ title(title_str)
%
% Solutions at the finest grid
%
dense=h5read(file,'/dense');
x=h5read(file,'/solutions/xg');
y=h5read(file,'/solutions/yg');
sol_anal=h5read(file,'/solutions/anal');
%sol_direct=h5read(file,'/solutions/direct');
sol_calc=h5read(file,'/solutions/calc');
figure
% $$$ surf(x,y,sol_calc'-sol_anal')
pcolor(x,y,sol_calc'-sol_anal')
shading interp
colorbar
xlabel('X'); ylabel('Y'); zlabel('Error');
title(title_str)
figure
subplot(211)
[yy,iy] = max(abs(sol_anal),[],2);
[xx,ix] = max(yy);
iy0=iy(ix);
str=sprintf('Solution at y = %.4f', y(iy0));
plot(x, sol_anal(:,iy0),x, sol_calc(:,iy0),'o')
xlabel('x'); ylabel(str);
grid on
legend('Analytic Solution', 'MG Solution')
title(title_str)
subplot(212)
[xx,ix] = max(abs(sol_anal));
[yy,iy] = max(xx);
ix0=ix(iy);
str=sprintf('Solution at x = %.4f', x(ix0));
plot(y, sol_anal(ix0,:),y, sol_calc(ix0,:),'o')
xlabel('y'); ylabel(str);
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