Skip to content
Snippets Groups Projects
Commit 9530ffa4 authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann
Browse files

matlab scripts for post processing

parent 5a5dc3cd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
matlab/default_plots_options.m 0 → 100644
+ 9
− 0
View file @ 9530ffa4
%% Default values for plots
set(0,'defaulttextInterpreter','latex')
set(0, 'defaultAxesTickLabelInterpreter','latex');
set(0, 'defaultLegendInterpreter','latex');
set(0,'defaultAxesFontSize',16)
set(0, 'DefaultLineLineWidth', 1.5);
line_colors = [[0, 0.4470, 0.7410];[0.8500, 0.3250, 0.0980];[0.9290, 0.6940, 0.1250];...
[0.4940, 0.1840, 0.5560];[0.4660, 0.6740, 0.1880];[0.3010, 0.7450, 0.9330];[0.6350, 0.0780, 0.1840]];
line_styles = {'-','-.','--',':'};
\ No newline at end of file
matlab/helaz_analysis.m 0 → 100644
+ 112
− 0
View file @ 9530ffa4
%% HeLaZ data
filename = 'results_00.h5';
default_plots_options % Script to set up default plot variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load the data
moment = 'Ni00';
kr = h5read(filename,['/data/var2d/' moment '/coordkr']);
kz = h5read(filename,['/data/var2d/' moment '/coordkz']);
timeNi = h5read(filename,'/data/var2d/time');
Nipj = zeros(numel(timeNi),numel(kr),numel(kz));
for it = 1:numel(timeNi)
tmp = h5read(filename,['/data/var2d/', moment,'/', num2str(it,'%06d')]);
Nipj(it,:,:) = tmp.real + 1i * tmp.imaginary;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot growth rate vs t
gammas = zeros(numel(kr),numel(kz));
shifts = zeros(numel(kr),numel(kz));
% Linear fit of log(Napj)
x1 = timeNi;
itmin = ceil(0.5 * numel(timeNi)); %Take the second half of the time evolution
for ikr = 1:numel(kr)
for ikz = 1:numel(kz)
fit = polyfit(x1(itmin:end),log(abs(Nipj(itmin:end,ikr,ikz))),1);
gammas(ikr,ikz) = fit(1);
shifts(ikr,ikz) = fit(2);
end
end
FIGNAME = 'gamma_t';
fig = figure;
for ikr = 1:numel(kr)
linename = ['$k_r = ',num2str(kr(ikr)),'$'];
plot(kz,gammas(ikr,:),'DisplayName',linename);
end
TITLE = [];
TITLE = [TITLE,'$\eta_n=',num2str(1.0/MODEL.eta_n),'$, '];
TITLE = [TITLE,'$\eta_B=',num2str(MODEL.eta_B),'$, '];
TITLE = [TITLE, '$\nu=',num2str(MODEL.nu),'$, '];
%TITLE = [TITLE, '$k_z=',num2str(GRID.kz),'$'];
title(TITLE);
grid on
legend('show')
xlabel('$t$')
ylabel(['$|',moment,'|$'])
%% Saving fig
if SAVEFIG
FIGDIR = ['../results/', SIMID,'/'];
if ~exist(FIGDIR, 'dir')
mkdir(FIGDIR)
end
FIGNAME = [FIGNAME,'_Pe_',num2str(GRID.pmaxe),'_Je_',num2str(GRID.jmaxe),...
'_Pi_',num2str(GRID.pmaxi),'_Ji_',num2str(GRID.jmaxi),...
'_etan_',num2str(MODEL.eta_n),'_etaB_',num2str(MODEL.eta_B),'_nu_',num2str(MODEL.nu)];
FIGNAME = [FIGDIR, FIGNAME,'.fig'];
savefig(fig,FIGNAME);
disp(['Figure saved @ : ',FIGNAME])
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot Ni00 evolution
FIGNAME = 'Ni00_t';
fig = figure;
%HeLaZ results
x1 = timeNi;
il = 1;
for ikr = 1:numel(kr)
ic = 1;
for ikz = 1:2:numel(kz)
linename = ['$k_r = ',num2str(kr(ikr)),', k_z = ', num2str(kz(ikz)),'$'];
y1 = abs(Nipj(:,ikr,ikz));
semilogy(x1,y1,...
'DisplayName',linename,'Color', line_colors(ic,:), 'LineStyle', '--')
hold on
semilogy(x1(itmin:end),exp(gammas(ikr,ikz)*x1(itmin:end) + shifts(ikr,ikz)),...
'DisplayName',[linename,' fit'],'Color', line_colors(ic,:), 'LineStyle', '-.')
ic = ic + 1;
end
il = il + 1;
end
TITLE = [];
TITLE = [TITLE,'$\eta_n=',num2str(1.0/MODEL.eta_n),'$, '];
TITLE = [TITLE,'$\eta_B=',num2str(MODEL.eta_B),'$, '];
TITLE = [TITLE, '$\nu=',num2str(MODEL.nu),'$, '];
%TITLE = [TITLE, '$k_z=',num2str(GRID.kz),'$'];
title(TITLE);
grid on
legend('show')
xlabel('$t$')
ylabel(['$|',moment,'|$'])
%% Saving fig
if SAVEFIG
FIGDIR = ['../results/', SIMID,'/'];
if ~exist(FIGDIR, 'dir')
mkdir(FIGDIR)
end
FIGNAME = [FIGNAME,'_Pe_',num2str(GRID.pmaxe),'_Je_',num2str(GRID.jmaxe),...
'_Pi_',num2str(GRID.pmaxi),'_Ji_',num2str(GRID.jmaxi),...
'_etan_',num2str(MODEL.eta_n),'_etaB_',num2str(MODEL.eta_B),'_nu_',num2str(MODEL.nu)];
FIGNAME = [FIGDIR, FIGNAME,'.fig'];
savefig(fig,FIGNAME);
disp(['Figure saved @ : ',FIGNAME])
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
matlab/write_fort90.m 0 → 100644
+ 58
− 0
View file @ 9530ffa4
function [INPUT] = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC)
% Write the input script "fort.90" with desired parameters
INPUT = '../wk/fort.90';
fid = fopen(INPUT,'wt');
fprintf(fid,'&BASIC\n');
fprintf(fid,[' nrun=', num2str(BASIC.nrun),'\n']);
fprintf(fid,[' dt=', num2str(BASIC.dt),'\n']);
fprintf(fid,[' tmax=', num2str(BASIC.tmax),' ! time normalized to 1/omega_pe\n']);
fprintf(fid,'/\n');
fprintf(fid,'&GRID\n');
fprintf(fid,[' pmaxe =', num2str(GRID.pmaxe),' ! number of Hermite moments \n']);
fprintf(fid,[' jmaxe = ', num2str(GRID.jmaxe),' ! number of Hermite moments \n']);
fprintf(fid,[' pmaxi = ', num2str(GRID.pmaxi),' ! number of Hermite moments \n']);
fprintf(fid,[' jmaxi = ', num2str(GRID.jmaxi),' ! number of Hermite moments\n']);
fprintf(fid,[' nkr = ', num2str(GRID.nkr),'\n']);
fprintf(fid,[' krmin = ', num2str(GRID.krmin),'\n']);
fprintf(fid,[' krmax = ', num2str(GRID.krmax),' ! Normalized to cs0/Omega_i\n']);
fprintf(fid,[' nkz = ', num2str(GRID.nkz),'\n']);
fprintf(fid,[' kzmin = ', num2str(GRID.kzmin),'\n']);
fprintf(fid,[' kzmax = ', num2str(GRID.kzmax),' ! Normalized to cs0/Omega_i\n']);
fprintf(fid,'/\n');
fprintf(fid,'&OUTPUT_PAR\n');
fprintf(fid,[' nsave_0d = ', num2str(OUTPUTS.nsave_0d),'\n']);
fprintf(fid,[' nsave_1d = ', num2str(OUTPUTS.nsave_1d),'\n']);
fprintf(fid,[' nsave_2d = ', num2str(OUTPUTS.nsave_2d),'\n']);
fprintf(fid,[' nsave_5d = ', num2str(OUTPUTS.nsave_5d),'\n']);
fprintf(fid,[' write_Ni00 = ', OUTPUTS.write_Ni00,'\n']);
fprintf(fid,[' write_moments = ', OUTPUTS.write_moments,'\n']);
fprintf(fid,[' write_phi = ', OUTPUTS.write_phi,'\n']);
fprintf(fid,[' write_doubleprecision = ', OUTPUTS.write_doubleprecision,'\n']);
fprintf(fid,[' resfile0 = ', OUTPUTS.resfile0,'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&MODEL_PAR\n');
fprintf(fid,' ! Collisionality\n');
fprintf(fid,[' nu = ', num2str(MODEL.nu),' ! Normalized collision frequency normalized to plasma frequency\n']);
fprintf(fid,[' tau_e = ', num2str(MODEL.tau_e),' ! T_e/T_e\n']);
fprintf(fid,[' tau_i = ', num2str(MODEL.tau_i),' ! T_i/T_e temperature ratio\n']);
fprintf(fid,[' sigma_e = ', num2str(MODEL.sigma_e),' ! sqrt(m_e/m_i) mass ratio\n']);
fprintf(fid,[' sigma_i = ', num2str(MODEL.sigma_i),' ! sqrt(m_i/m_i)\n']);
fprintf(fid,[' q_e = ', num2str(MODEL.q_e),' ! Electrons charge\n']);
fprintf(fid,[' q_i = ', num2str(MODEL.q_i),' ! Ions charge\n']);
fprintf(fid,[' eta_n = ', num2str(MODEL.eta_n),' ! L_perp / L_n Density gradient\n']);
fprintf(fid,[' eta_T = ', num2str(MODEL.eta_T),' ! L_perp / L_T Temperature gradient\n']);
fprintf(fid,[' eta_B = ', num2str(MODEL.eta_B),' ! L_perp / L_B Magnetic gradient and curvature\n']);
fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),' ! Debye length\n']);
fprintf(fid,'/\n');
fprintf(fid,'&INITIAL_CON\n');
fprintf(fid,' ! Background value\n');
fprintf(fid,[' initback_moments=', num2str(INITIAL.initback_moments),' ! x 1e-3\n']);
fprintf(fid,' ! Noise amplitude\n');
fprintf(fid,[' initnoise_moments=', num2str(INITIAL.initnoise_moments),'\n']);
fprintf(fid,[' iseed=', num2str(INITIAL.iseed),'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&TIME_INTEGRATION_PAR\n');
fprintf(fid,[' numerical_scheme=', TIME_INTEGRATION.numerical_scheme,'\n']);
fprintf(fid,'/');
fclose(fid);
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment