Homec4science
Diffusion SPClibs (feature/fft-split) Browse

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

Last Change
Raw File

test_relax.m
View Options

%
% @file test_relax.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>
%
file='test_relax.h5';
%
nx=h5readatt(file,'/','NX');
alpha=h5readatt(file,'/','ALPHA');
nidbas=h5readatt(file,'/','NIDBAS');
relax=h5readatt(file,'/','RELAX');
omega=h5readatt(file,'/','OMEGA');
if alpha == 0
kx=h5readatt(file,'/','KX');
title_str=sprintf('N=%d, NIDBAS=%d, KX=%d, relax=%s, omega=%.3f', ...
nx,nidbas,kx,relax, omega);
else
modem=h5readatt(file,'/','MODEM');
modep=h5readatt(file,'/','MODEP');
title_str=sprintf('N=%d, NIDBAS=%d, modem=%d, modep=%d, relax=%s, omega=%.3f', ...
nx,nidbas,modem,modep,relax, omega);
end
%
% Solutions at the finest grid
%
x=h5read(file,'/solutions/xg');
sol_direct=h5read(file,'/solutions/direct');
sol_anal=h5read(file,'/solutions/anal');
sol_calc=h5read(file,'/solutions/calc');
figure
subplot(211)
plot(x, sol_anal, x, sol_calc, 'o')
legend('Analytic', 'Calculated')
xlabel('X')
grid on
title(title_str);
%
% Relaxations
%
errdisc=h5read(file,'/relaxation/errdisc');
resid=h5read(file,'/relaxation/resid');
its=0:1:size(errdisc)-1;
subplot(212)
semilogy(its,errdisc,its,resid)
legend('Discretisation Error', 'Residual Norm')
xlabel('Iterations')
grid on
title(title_str)
%
% FE Matrix
%
dset = '/MATA/';
A = gb_mat(file, dset);
D = diag(diag(A),0);
n=rank(A);
k=1:1:n;
if relax(1:2) == 'ja'
%
% Compute eigenvalues of Rj = D^(-1)*A
%
[V, l] = eig(A,D);
[lambda, iss] = sort(diag(l));
V = V(1:end,iss);
%
% Spectral radius of Jacobi iteration matrix
% R(omega) = max |1-omega*lambda|
%
om=0:0.01:1;
for i=1:size(om,2)
rho(i) = max(abs(1-om(i).*lambda));
end
fprintf(1, 'Spectral Radius of relaxation matrix = %g\n', max(abs(1-omega*lambda)))
figure
subplot(211)
plot(k, 1-omega*lambda, 'o-')
xlabel('mode k'); ylabel('Eigen value of inv(D)*A')
grid on
title(title_str)
subplot(212)
plot(om, rho)
omega_c = 2.0/max(lambda);
str = sprintf('Critical omega = %.3f', omega_c)
title(str)
xlabel('\omega'); ylabel('Spectral Radius')
grid on
elseif relax(1:2) == 'gs'
%
% Spectral radius of GS Iteration Matrix
% Rg = -(D+L)^(-1) * U
%
B = tril(A,0); % D+L
% $$$ [V, l] = eig(-triu(A,1),B); lambda=diag(l);
% $$$ [V, l] = eig(B,A); lambda = 1 - 1./diag(l);
[V, l] = eig(A,B); lambda = 1 - diag(l);
[lambda, iss] = sort(lambda, 'descend');
V = V(1:end,iss);
fprintf(1, 'Spectral Radius of relaxation matrix = %g\n', max(abs(lambda)))
figure
subplot(211)
plot(real(lambda), imag(lambda), 'o')
xlabel('Real of eigenvalues'); ylabel(['Imag of ' ...
'eigenvalues'])
axis equal
title(title_str)
grid on
subplot(212)
plot(k, abs(lambda), 'o-')
xlabel('Mode'); ylabel('eigenvalues')
grid on
end
%
% Plot eigenvectors
neig=size(lambda,1);
for i=1:neig
k=mod(i-1,4*5)+1;
if k==1
figure
title(title_str)
end
subplot(4,5,k)
if relax(1:2) == 'ja'
str = sprintf('Mode = %d, ||R|| = %.4f', i, 1-omega* ...
lambda(i));
else
str = sprintf('Mode = %d, ||R|| = %.4f', i, ...
lambda(i));
end
plot(V(:,i)); grid on
title(str)
end
c4science · Help