Homec4science
Diffusion SPClibs (feature/fft-split) Browse

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

Last Change
Raw File

pde2d_sym_pardiso_dft.m
View Options

%
% @file pde2d_sym_pardiso_dft.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>
%
mat_disp=1;
file='pde2d_sym_pardiso_dft.h5';
if verLessThan('matlab', '7.9');
n = hdf5read(file,'/MAT1/RANK');
nnz = hdf5read(file,'/MAT1/NNZ');
nlsym = hdf5read(file,'/MAT1/NLSYM');
NR=hdf5read(file,'/NX'); NTH=hdf5read(file,'/NY');
NIDBAS1=hdf5read(file,'/NIDBAS1');
NIDBAS2=hdf5read(file,'/NIDBAS2');
MBESS=hdf5read(file,'/MBESS');
EPSI=hdf5read(file,'/EPSI');
KMIN=hdf5read(file,'/KMIN');
KMAX=hdf5read(file,'/KMAX');
else
n = hdf5read(file,'/MAT1/', 'RANK');
nnz = hdf5read(file,'/MAT1/', 'NNZ');
nlsym = hdf5read(file,'/MAT1/', 'NLSYM');
NR=hdf5read(file,'/', 'NX'); NTH=hdf5read(file,'/', 'NY');
NIDBAS1=hdf5read(file,'/', 'NIDBAS1');
NIDBAS2=hdf5read(file,'/', 'NIDBAS2');
MBESS=hdf5read(file,'/', 'MBESS');
EPSI=hdf5read(file,'/','EPSI');
KMIN=hdf5read(file,'/','KMIN');
KMAX=hdf5read(file,'/','KMAX');
end
LABEL=sprintf('nr = %d, ntheta = %d, nidbas = (%d,%d), mbess = %d, epsi = %3.2f', NR, NTH, ...
NIDBAS1, NIDBAS2, MBESS, EPSI);
DK = KMAX-KMIN+1;
attr=hdf5read(file,'/xgrid/','title'); x_ann=attr.Data;
attr=hdf5read(file,'/ygrid/','title'); y_ann=attr.Data;
attr=hdf5read(file,'/sol/','title'); sol_ann=attr.Data;
attr=hdf5read(file,'/solana/','title'); solexact_ann=attr.Data;
attr=hdf5read(file,'/errors/','title');err_ann=attr.Data;
if mat_disp == 1
cols=hdf5read(file, '/MAT1/cols');
irow=hdf5read(file, '/MAT1/irow');
val=h5Complex(file, '/MAT1/val');
perm=hdf5read(file, '/MAT1/perm');
rows = zeros(nnz,1);
cols=double(cols);
irow=double(irow);
perm=double(perm);
for i=1:n
s = irow(i);
e = irow(i+1)-1;
rows(s:e) = i;
end
valr=real(val); vali=imag(val);
mat = sparse(rows,cols,valr);
figure
subplot(121)
spy(mat,8)
title('Original Matrix structure')
subplot(122)
spy(mat(perm,perm),8)
title('Permuted Matrix structure')
end
r=hdf5read(file,'/xgrid');
t=hdf5read(file,'/ygrid');
sol=hdf5read(file,'/sol')';
solexact=hdf5read(file,'/solana')';
err=hdf5read(file,'/errors')';
solr=hdf5read(file,'/derivx')';
solt=hdf5read(file,'/derivy')';
[R,T]=meshgrid(r,t);
x = R.*cos(T); y= R.*sin(T);
solx = cos(T).*solr - sin(T)./R.*solt;
soly = sin(T).*solr + cos(T)./R.*solt;
figure
subplot(221)
pcolor(r,t,sol);
shading interp
hold on, quiver(r,t,solr,solt)
xlabel(x_ann); ylabel(y_ann)
title(LABEL)
colorbar
subplot(222)
pcolor(x,y,sol)
shading interp
hold on, quiver(x,y,solx,soly)
hold off, axis image
xlabel('X'); ylabel('Y')
title('X-Y plane')
colorbar
subplot(223)
surfc(x,y,sol)
xlabel('X'); ylabel('Y');
title(sol_ann)
subplot(224)
surfc(x,y,err)
xlabel('X'); ylabel('Y');
title(err_ann)
figure
ft_sol=h5Complex(file,'/FT_SOL');
ft_sol=reshape(ft_sol,DK,[]);
m=[KMIN:KMAX]; sp=1:NR+NIDBAS1;
subplot(121)
stem3(sp,m, real(ft_sol), 'filled')
shading interp
xlabel('Radial spline number'); ylabel('m')
title('Real(\phi)')
subplot(122)
stem3(sp,m, imag(ft_sol),'filled')
shading interp
xlabel('Radial spline number'); ylabel('m')
title('Imag(\phi)')
figure
ft_rhs=h5Complex(file,'/FT_RHS');
ft_rhs=reshape(ft_rhs,DK,[]);
m=[KMIN:KMAX]; sp=1:NR+NIDBAS1;
subplot(121)
stem3(sp,m, real(ft_rhs), 'filled')
shading interp
xlabel('Radial spline number'); ylabel('m')
title('Real(\rho)')
subplot(122)
stem3(sp,m, imag(ft_rhs),'filled')
shading interp
xlabel('Radial spline number'); ylabel('m')
title('Imag(\rho)')
figure
energy_k = h5Complex(file,'/ENERGY_K');
subplot(211)
stem(m, real(energy_k));
xlabel('m'); ylabel('Real(\phi)'); title('Spectral energy')
subplot(212)
stem(m, imag(energy_k));
xlabel('m'); ylabel('Imag(\phi)'); title('Spectral energy')
c4science · Help