Homec4science
Diffusion yearII_reports (master) Browse

rSTICAZZI/physnum/rap6/Conv_nontrivial.m
fe2a2725821dmaster

Last Change
Raw File

Conv_nontrivial.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 = 'exercice6'; % Nom de l'executable (NB: ajouter .exe sous Windows)
input = 'configuration.in'; % Nom du fichier d'entree de base
nsimul = 501; % Nombre de simulations a faire
coeff = 4.0/5;
N_1 = linspace(10, nsimul+9, nsimul);
N_2 = floor(coeff * N_1);
outfile = sprintf('outs/conv_coeff=%fnontrival.dat', coeff);
outE = cell(1, nsimul);
outdir = 'm_outs/';
p = 1;
chooseN1 = [30, 50, 80, 120, 150, 180, 200, 250, 350, 400, 480];
%% Simulations %%
%%%%%%%%%%%%%%%%%
output = cell(1, nsimul); % Tableau de cellules contenant le nom des fichiers de sortie
interval = cell(1, nsimul);
for i = 1:nsimul
output{i} = [outdir, 'N=', num2str(N_1(i)), 'conv_nontrivial.out'];
outE{i} = ['outs/field/', num2str(N_1(i)), 'Efield.dat'];
interval{i} = sprintf('{{0.02, %d}, {0.12, %d}}', N_1(i), N_2(i));
% Execution du programme en lui envoyant la valeur a scanner en argument
cmd = sprintf('%s%s %s \"interval=%s\" output=%s trapezium=%f trivial=false', repertoire, executable, input, interval{i}, output{i}, p);
disp(cmd);
system(cmd);
end
%% Analyse %%
%%%%%%%%%%%%%
% Parcours des resultats de toutes les simulations
a0 = 1e4;
R = 0.12;
b = 0.02;
phi_err = zeros(1,nsimul);
chooseIndex = 1;
for i = 1:nsimul % Parcours des resultats de toutes les simulations
data = load(output{i});
[rmid, E] = electric_field(data(:,1), data(:,3));
if chooseIndex < length(chooseN1) && N_1(i) == chooseN1(chooseIndex)
% write electric field on files
eout = fopen(outE{i}, 'w');
fprintf(eout, 'rmid E r phi\n');
for j = 1:length(rmid)
fprintf(eout, '%.15g %.15g %.15g %.15g %.15g\n', rmid(j), E(j), E(j) * rmid(j) / epsilon(rmid(j), b, R), data(j,1), data(j,3));
end
fclose(eout);
chooseIndex = chooseIndex + 1;
end
phi_err(i) = data(N_1(i)+1,3);
fprintf('N = %d, phi(b) = %.15g\n', N_1(i) + N_2(i), data(N_1(i)+1,3));
end
phi_conv = phi_err(nsimul);
phi_err = abs(phi_err - phi_conv);
fprintf('Convergence value: %f\n', phi_conv);
%% Linear regression %%
%%%%%%%%%%%%%%%%%%%%%%%
m = polyfit(log(N_1), log(phi_err), 1);
linreg = fprintf('Loglog fit: %.8g * x + %.8g', m(1), m(2));
%% File writing %%
%%%%%%%%%%%%%
% print data on file
out = fopen(outfile, 'w');
fprintf(out, 'N1 N2 phi_err\n');
for i = 1:nsimul
fprintf(out, '%d %d %.15g\n', N_1(i), N_2(i), phi_err(i));
end
fclose(out);
%% File writing %%
%%%%%%%%%%%%%
% print data on file
out = fopen(outfile, 'w');
fprintf(out, 'N1 N2 phi_err\n');
for i = 1:nsimul
fprintf(out, '%d %d %.15g\n', N_1(i), N_2(i), phi_err(i));
end
fclose(out);
%% Figures %%
%%%%%%%%%%%%%
% make figure with data
figure
loglog(N_1, phi_err,'k+')
xlabel('N_1')
ylabel('phi_err')
grid on
%% Helper %%
%%%%%%%%%%%%
function [rmid, E] = electric_field(r, phi)
N = length(r) - 1;
rmid = zeros(N, 1);
E = zeros(N, 1);
for i=1:N
h = r(i+1) - r(i);
rmid(i) = r(i) + h;
E(i) = - (phi(i+1) - phi(i))/ h;
end
end
function eps = epsilon(r, b, R)
if (r > b && abs(r-b) > 1e-13)
eps = 8 - 6 * (r-b)/(R-b);
else
eps = 1;
end
end
c4science · Help