Homec4science
Diffusion yearII_reports (master) Browse

rSTICAZZI/physnum/rap8/convergence.m
e0dc37e4725amaster

Last Change
Raw File

convergence.m
View Options

% Ce script Matlab automatise la production de resultats
% lorsqu'on doit faire une serie de simulations en
% variant un des parametres d'entree.
%
% Il utilise les arguments du programme (voir ConfigFile.h)
% pour remplacer la valeur d'un parametre du fichier d'input
% par la valeur scannee.
%
%% Parametres %%
%%%%%%%%%%%%%%%%
repertoire = './'; % Chemin d'acces au code compile (NB: enlever le ./ sous Windows)
executable = 'exercice8'; % Nom de l'executable (NB: ajouter .exe sous Windows)
input = 'configuration.in'; % Nom du fichier d'entree de base
tfin = 5000;
xl = -200;
xr = 200;
bounds = sprintf('(%.15g, %.15g)', xl, xr);
omega = 0.003;
delta = 0;
x0 = 0;
sigma = 0.06;
n = 14;
N = 300;
hbar = 1.0;
m = 1.0;
n_sym = 50;
sampling = 200;
%dt = logspace(-1, 0, n_sym);
dt = zeros(n_sym, 1);
for i=1:n_sym
dt(i) = tfin / (1000 * i + 500);
end
dtlast = 1e-2;
dt(n_sym+1) = dtlast;
output = cell(1, n_sym+1);
x = gen_space(xl, xr, N);
xmid = zeros(n_sym+1, 1);
pmid = zeros(n_sym+1, 1);
dxmid = zeros(n_sym+1, 1);
dpmid = zeros(n_sym+1, 1);
p0 = hbar * 2 * pi * n / (xr - xl);
A = p0 / (omega * m);
%% Simulations %%
%%%%%%%%%%%%%%%%%
for i=1:(n_sym+1)
output{i} = ['m_outs/conv_', num2str(dt(i))];
cmd = sprintf('%s%s %s N=%d tfin=%.15g dt=%.15g "bounds=%s" omega=%.15g x0=%.15g sigma=%.15g n=%.15g delta=%.15g output=%s sampling=%d outmode=%d', ...
repertoire, executable, input,...
N, tfin, dt(i), ...
bounds, omega, x0, sigma, n, delta, ...
output{i}, sampling, 1);
disp(cmd);
system(cmd);
state = load(sprintf('%s_state.out', output{i}));
xmid(i) = state(end, 5);
dxmid(i) = state(end, 7);
pmid(i) = state(end, 8);
dpmid(i) = state(end, 10);
end
%% Plot %%
%xref = A * sin(omega * tfin);
%pref = p0 * cos(omega * tfin);
xref = xmid(end);
pref = pmid(end);
dxref = dxmid(end);
dpref = dpmid(end);
dt_ = dt(1:end-1);
diffx = abs(xmid(1:end-1) - xref);
figure
loglog(dt_, diffx, 'k+')
%plot(dt_.^2, diffx, 'k+')
%plot(dt.^2, xmid, 'k+')
grid on
xlabel('dt')
ylabel('<x>')
fitx = polyfit(log(dt_), log(diffx), 1);
diffdx = abs(dxmid(1:end-1) - dxref);
figure
loglog(dt_, diffdx, 'k+')
%plot(dt_.^2, diffx, 'k+')
%plot(dt.^2, xmid, 'k+')
grid on
xlabel('dt')
ylabel('<dx>')
fitdx = polyfit(log(dt_), log(diffdx), 1);
diffp = abs(pmid(1:end-1) - pref);
figure
loglog(dt_, diffp, 'k+')
%plot(dt.^2, pmid, 'k+')
grid on
xlabel('dt')
ylabel('<p>')
fitp = polyfit(log(dt_), log(diffp), 1);
diffdp = abs(dpmid(1:end-1) - dpref);
figure
loglog(dt_, diffdp, 'k+')
%plot(dt_.^2, diffx, 'k+')
%plot(dt.^2, xmid, 'k+')
grid on
xlabel('dt')
ylabel('<dp>')
fitdp = polyfit(log(dt_), log(diffdp), 1);
%figure
%grid on
%xlabel('dt')
%ylabel('<dx>')
%plot(dt, abs(dxmid - dxref), 'k+')
%figure
%grid on
%xlabel('dt')
%ylabel('<dp>')
%plot(dt, dpmid, 'k+')
%% Print data %%
outfile = 'outs/conv.dat';
% print data on file
out = fopen(outfile, 'w');
fprintf(out, 'dt x p dx dp\n');
for i = 1:n_sym
fprintf(out, '%.15g %.15g %.15g %.15g %.15g\n', dt_(i), diffx(i), diffp(i), diffdx(i), diffdp(i));
end
fclose(out);
%% Utils %%
function x = gen_space(xl, xr, N)
x = zeros(N+1, 1);
dx = (xr - xl) / N;
for k = 0:N
x(k+1) = k * dx + xl;
end
end
c4science · Help