Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/matlab/jac_opt.m
93a1e9712d04feature/fft-split

Last Change
Raw File

jac_opt.m
View Options

%
% @file jac_opt.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>
%
%
% Find optimal damping factor \omega for Jacobi relaxations
%
clear all
alpha=0.5;
theta2=pi/2:0.01:pi;
theta1=0:0.01:pi;
[x,y]=meshgrid(theta1,theta1);
sint12=sin(theta1./2).^2;
sint22=sin(theta2./2).^2;
[ksi,eta]=meshgrid(sint12,sint22);
nth1=length(theta1);
nth2=length(theta2);
omega=0.1:0.01:1;
n=length(omega);
for i=1:n
c=2*omega(i)/(1+alpha^2);
G=abs(1 - c.*(ksi+(alpha^2).*eta ));
[gmax,imax]=max(G);
[mu(i),jmax]=max(gmax);
eta_opt(i)=eta(imax(jmax),jmax);
ksi_opt(i)=ksi(imax(jmax),jmax);
end
[mu_min,i_min]=min(mu);
omega_opt=omega(i_min);
str_title=sprintf('omega = %.3f, mu = %.3f, alpha = %.2f', omega_opt, ...
mu_min, alpha);
figure
plot(omega,mu,'o-', omega, ksi_opt, '*-', omega, eta_opt, '^-');
legend('\mu', '\xi_{opt}', '\eta_{opt}')
xlabel('\omega'); ylabel('\mu')
grid on
title(str_title);
% $$$
% $$$ c=2*omega_opt/(1+alpha^2);
% $$$ G=1-c.*(ksi+(alpha^2).*eta);
% $$$ figure
% $$$ hold off
% $$$ plot(theta1, G(1,:), 'r', 'LineWidth', 2)
% $$$ hold on
% $$$ plot(theta1, G(nth2,:), 'g', 'LineWidth', 2)
% $$$ for jj=1:20:nth2
% $$$ plot(theta1, G(jj,:), 'b')
% $$$ end
% $$$ xlabel('\theta_1'); ylabel('Amplification Factor for Jacobi')
% $$$ title(str_title)
% $$$ grid on
c4science · Help