Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/matlab/tcdsmat.m
3324689e437cfeature/refactoring

Last Change
Raw File

tcdsmat.m
View Options

%
% @file tcdsmat.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>
%
mat='/MAT1';
gbmat;
clear S gb_mat;
file='tcdsmat.h5';
nx=hdf5read(file,'/NX');
ny=hdf5read(file,'/NY');
dists=hdf5read(file,strcat(mat,'/dists'));
vals=hdf5read(file,strcat(mat,'/vals'));
rowv=hdf5read(file,strcat(mat,'/rowv'));
colh=hdf5read(file,strcat(mat,'/colh'));
n=hdf5read(file,strcat(mat,'/RANK'));
nd=hdf5read(file,strcat(mat,'/NDIAGS'));
err=zeros(n,nd);
% Diagonal balancing of matrix
dbal = 1./sqrt(diag(a));
a = diag(dbal)*a*diag(dbal);
% Check CDS mat except row ny and column ny
for k=1:nd
d=dists(k);
i1=max(1,1-d); i2=min(n,n-d);
fprintf(1,'%8d %8d %8d\n',d,i1,i2);
for i=i1:i2
j=i+d;
if (i~=ny && j~=ny)
err(i,k) = a(i,j)-vals(i,k);
end
end
end
fprintf(1,'min/max of err: %8.4e, %8.4e\n',min(min(err)), max(max(err)));
% Check row ny and j .ne. ny
i=ny;
bw0=size(rowv,1);
for k=1:nd
d=dists(k);
j=i+d;
if ((j >= ny+1) && (j <= bw0))
err(i,k)=a(i,j)-rowv(j);
end
end
fprintf(1,'min/max of err: %8.4e, %8.4e\n',min(min(err)), max(max(err)));
% Check column ny
j=ny;
for k=1:nd
d=dists(k);
i=j-d;
if ((i >= ny+1) && (i <= bw0))
err(i,k)=a(i,j)-colh(i);
end
end
fprintf(1,'min/max of err: %8.4e, %8.4e\n',min(min(err)), max(max(err)));
% Check RHS
rhs=hdf5read(file,'/RHS');
fprintf('Err in RHS: %8.3e\n', max(max(abs(rhs-rhs0))))
% Check SOL
sol=hdf5read(file,'/SOL');
err= sol-sol0;
fprintf('Err SOL: %8.3e\n', max(max(abs(err))));
c4science · Help