Homec4science
Diffusion SPClibs (feature/fft-split) Browse

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

Last Change
Raw File

test_mg.m
View Options

%
% @file test_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>
%
file='test_mg.h5';
%
nx=h5readatt(file,'/','NX');
nidbas=h5readatt(file,'/','NIDBAS');
relax=h5readatt(file,'/','RELAX');
levels=h5readatt(file,'/','LEVELS');
nu1=h5readatt(file,'/','NU1');
nu2=h5readatt(file,'/','NU2');
alpha=h5readatt(file,'/','ALPHA');
omega=h5readatt(file,'/','OMEGA');
if alpha == 0
kx=h5readatt(file,'/','KMODE');
title_str=sprintf('N=%d, NIDBAS=%d, KX=%d, relax=%s, omega=%.3f, levels = %d, nu1 = %d, nu2 = %d', ...
nx,nidbas,kx,relax, omega, levels, nu1, nu2);
else
modem=h5readatt(file,'/','MODEM');
modep=h5readatt(file,'/','MODEP');
title_str=sprintf('N=%d, NIDBAS=%d, modem=%d, modep=%d, relax=%s, omega=%.3f, levels = %d, nu1 = %d, nu2 = %d', ...
nx,nidbas,modem,modep,relax, omega, levels, nu1, nu2);
end
%
% Read matrices at coarset grid
%
for lev=2:levels
%
% FE mat at fine grid
mglevel=sprintf('/mglevels/level.%.2d', lev-1);
dset=strcat(mglevel,'/mata');
mata_f = gb_mat(file, dset);
%
% FE mat at coarse grid
mglevel=sprintf('/mglevels/level.%.2d', lev);
dset=strcat(mglevel,'/mata');
mata_c = gb_mat(file, dset);
%
% Prolong mat
dset=strcat(mglevel,'/matp');
matp=h5read(file,dset);
%
% Check
fprintf(1,'Level %d: ||A_coarse - P''*A_fine*P|| = %g\n', lev, norm(matp'*mata_f*matp ...
- mata_c))
end
%
% Iterations
dset='/Iterations/';
err=h5read(file, strcat(dset,'errors'));
disc_err=h5read(file, strcat(dset,'disc_errors'));
resid=h5read(file, strcat(dset,'residues'));
its=0:1:size(err,1)-1;
figure
subplot(212)
semilogy(its,resid,'o-', its, disc_err,'h-')
legend('Residue', 'Error')
grid on
xlabel('Iterations'); ylabel('Norm od residue and error');
title(title_str);
%
% Plot grid values at the last iteration
xgrid=h5read(file, '/Iterations/xgrid');
u_calc=h5read(file, '/Iterations/u_calc');
u_exact=h5read(file, '/Iterations/u_exact');
u_direct=h5read(file, '/Iterations/u_direct');
subplot(211)
plot(xgrid, u_exact, xgrid,u_calc,'o')
legend('Analytic', 'Calculated')
xlabel('X');ylabel('Grid values of solution')
grid on
title(title_str);
%
%
% $$$ mglevel=sprintf('/mglevels/level.%.2d', 1);
% $$$ dset=strcat(mglevel,'/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)
% $$$ xlabel('\omega'); ylabel('Spectral Radius')
% $$$ grid on
% $$$
% $$$ for i=1:n
% $$$ k=mod(i-1,4*5)+1;
% $$$ if k==1
% $$$ figure
% $$$ title(title_str)
% $$$ end
% $$$ subplot(4,5,k)
% $$$ str = sprintf('Mode = %d, ||R|| = %.3f', i, 1-omega*lambda(i));
% $$$ plot(V(:,i)); grid on
% $$$ title(str)
% $$$ end
% $$$ elseif relax(1:2) == 'gs'
% $$$ %
% $$$ % Spectral radius of GS Iteration Matrix
% $$$ % Rg = (D-L)^(-1) * U
% $$$ %
% $$$ B = tril(A,0); % D-L
% $$$ lambda = eig(-triu(A,1),B);
% $$$ fprintf(1, 'Spectral Radius of relaxation matrix = %g\n', max(abs(lambda)))
% $$$ figure
% $$$ plot(lambda, 'o', 'MarkerSize', 6)
% $$$ axis equal
% $$$ grid on
% $$$ xlabel('Real of eigenvalues'); ylabel('Imag of eigenvalues')
% $$$ title(title_str)
% $$$ end
c4science · Help