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

Cleaned and refreshed

parent 3391be78
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
wk/benchmark_kperp_scan.m
+ 43
− 37
View file @ ef30a980
...@@ -12,19 +12,19 @@ OUTPUTS.write_phi = '.true.'; ...@@ -12,19 +12,19 @@ OUTPUTS.write_phi = '.true.';
OUTPUTS.write_doubleprecision = '.true.'; OUTPUTS.write_doubleprecision = '.true.';
OUTPUTS.resfile0 = '''results'''; OUTPUTS.resfile0 = '''results''';
%% Grid parameters %% Grid parameters
GRID.pmaxe = 15; GRID.pmaxe = 12;
GRID.jmaxe = 6; GRID.jmaxe = 7;
GRID.pmaxi = 15; GRID.pmaxi = 12;
GRID.jmaxi = 6; GRID.jmaxi = 7;
GRID.nkr = 1; GRID.nkr = 1;
GRID.krmin = 0.; GRID.krmin = 0.;
GRID.krmax = 0.; GRID.krmax = 0.;
GRID.nkz = 20; GRID.nkz = 20;
GRID.kzmin = 0.1; GRID.kzmin = 0.1;
GRID.kzmax = 1.5; GRID.kzmax = 2.0;
%% Model parameters %% Model parameters
MODEL.CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty) MODEL.CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
MODEL.nu = 0.1; % collisionality nu*L_perp/Cs0 MODEL.nu = 0.01; % collisionality nu*L_perp/Cs0
% temperature ratio T_a/T_e % temperature ratio T_a/T_e
MODEL.tau_e = 1.0; MODEL.tau_e = 1.0;
MODEL.tau_i = 1.0; MODEL.tau_i = 1.0;
...@@ -48,7 +48,18 @@ BASIC.tmax = 100.0; ...@@ -48,7 +48,18 @@ BASIC.tmax = 100.0;
INITIAL.initback_moments = 0.01; INITIAL.initback_moments = 0.01;
INITIAL.initnoise_moments = 0.; INITIAL.initnoise_moments = 0.;
INITIAL.iseed = 42; INITIAL.iseed = 42;
INITIAL.selfmat_file = ...
['''../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_',num2str(GRID.pmaxe),...
'_Jmaxe_',num2str(GRID.jmaxe),'_Pmaxi_',num2str(GRID.pmaxi),'_Jmaxi_',...
num2str(GRID.jmaxi),'_pamaxx_10.h5'''];
INITIAL.eimat_file = ...
['''../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_',num2str(GRID.pmaxe),...
'_Jmaxe_',num2str(GRID.jmaxe),'_Pmaxi_',num2str(GRID.pmaxi),'_Jmaxi_',...
num2str(GRID.jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
INITIAL.iemat_file = ...
['''../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_',num2str(GRID.pmaxe),...
'_Jmaxe_',num2str(GRID.jmaxe),'_Pmaxi_',num2str(GRID.pmaxi),'_Jmaxi_',...
num2str(GRID.jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Write input file %% Write input file
INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC); INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
...@@ -72,63 +83,58 @@ end ...@@ -72,63 +83,58 @@ end
%% Analysis and basic figures %% Analysis and basic figures
SAVEFIG = 1; SAVEFIG = 1;
filename = 'results_00.h5'; filename = 'results_00.h5';
default_plots_options
if MODEL.CO == -1; CONAME = 'FC';
elseif MODEL.CO == -2; CONAME = 'DC';
elseif MODEL.CO == 0; CONAME = 'LB'; end;
TITLE = [];
TITLE = [TITLE,'$\eta_n=',num2str(MODEL.eta_n),'$, '];
TITLE = [TITLE,'$\eta_B=',num2str(MODEL.eta_B),'$, '];
TITLE = [TITLE,'$\eta_T=',num2str(MODEL.eta_T),'$, '];
TITLE = [TITLE, '$\nu=',num2str(MODEL.nu),'$, '];
TITLE = [TITLE, '$(P,J)=(',num2str(GRID.pmaxe),',',num2str(GRID.jmaxe),')$'];
%% Growth rate analysis %% Growth rate analysis
gammas = zeros(numel(kr),numel(kz));
shifts = zeros(numel(kr),numel(kz));
moment = 'Ni00'; moment = 'Ni00';
kr = h5read(filename,['/data/var2d/' moment '/coordkr']); kr = h5read(filename,['/data/var2d/' moment '/coordkr']);
kz = h5read(filename,['/data/var2d/' moment '/coordkz']); kz = h5read(filename,['/data/var2d/' moment '/coordkz']);
timeNi = h5read(filename,'/data/var2d/time'); timeNi = squeeze(h5read(filename,'/data/var2d/time'));
Nipj = zeros(numel(timeNi),numel(kr),numel(kz)); Ni00 = zeros(numel(kr),numel(kz),numel(timeNi));
for it = 1:numel(timeNi) for it = 1:numel(timeNi)
tmp = h5read(filename,['/data/var2d/', moment,'/', num2str(it,'%06d')]); tmp = h5read(filename,['/data/var2d/', moment,'/', num2str(it,'%06d')]);
Nipj(it,:,:) = tmp.real + 1i * tmp.imaginary; Ni00(:,:,it) = tmp.real + 1i * tmp.imaginary;
end end
% Linear fit of log(Napj) % Amplitude ratio method
x1 = timeNi;
itmin = ceil(0.5 * numel(timeNi)); %Take the second half of the time evolution
ikr = 1; ikr = 1;
gammas = zeros(numel(kr),numel(kz));
Napj = ones(numel(timeNi),(GRID.pmaxe+1)*(GRID.jmaxe+1)+(GRID.pmaxi+1)*(GRID.jmaxi+1));
for ikz = 1:numel(kz) for ikz = 1:numel(kz)
fit = polyfit(x1(itmin:end),log(abs(Nipj(itmin:end,ikr,ikz))),1); gammas(ikr,ikz) = LinearFit_s(squeeze(timeNi), squeeze(abs(Ni00(ikr,ikz,:))));
gammas(ikr,ikz) = fit(1);
shifts(ikr,ikz) = fit(2);
end end
%% Plot %% Plot
fig = figure; fig = figure;
X = kz*sqrt(1+MODEL.tau_i); % Convert to Ricci 2006 Normalization
subplot(121) % growth rate
%HeLaZ results %HeLaZ results
X = kz*sqrt(1+MODEL.tau_i); Y1 = gammas(1,:);
Y = gammas(1,:); plot(X,Y1,'-','DisplayName','HeLaZ')
plot(X,Y,'-','DisplayName','HeLaZ')
hold on hold on
%MOLI results %MOLI results
X = kz*sqrt(1+MODEL.tau_i); Y2 = results.kperp.Maxgammas;
Y = results.kperp.Maxgammas; plot(X,Y2,'--','DisplayName','MOLI');
plot(X,Y,'--','DisplayName','MOLI');
title(TITLE); title(TITLE);
grid on grid on
legend('show') legend('show')
xlabel('$k_\perp * (1+\tau)^{1/2}$') xlabel('$k_\perp * (1+\tau)^{1/2}$')
ylabel('$\gamma L_\perp/c_{s} $') ylabel('$\gamma L_\perp/c_{s} $')
subplot(122) % Error
ERR = abs(gammas(1,:)' - results.kperp.Maxgammas);
plot(X,ERR,'-','DisplayName','MOLI relative error');
grid on
legend('show')
xlabel('$k_\perp * (1+\tau)^{1/2}$')
ylabel('$\epsilon_\gamma $')
if SAVEFIG; FIGNAME = 'g_vs_kz'; save_figure; end; if SAVEFIG; FIGNAME = 'g_vs_kz'; save_figure; end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
wk/benchmark_time_solver.m
+ 59
− 82
View file @ ef30a980
...@@ -21,11 +21,11 @@ GRID.nkr = 1; ...@@ -21,11 +21,11 @@ GRID.nkr = 1;
GRID.krmin = 0.; GRID.krmin = 0.;
GRID.krmax = 0.; GRID.krmax = 0.;
GRID.nkz = 1; GRID.nkz = 1;
GRID.kzmin = 1.0; GRID.kzmin = 0.67;
GRID.kzmax = 1.0; GRID.kzmax = 0.67;
%% Model parameters %% Model parameters
MODEL.CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty) MODEL.CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
MODEL.nu = 1.0; % collisionality nu*L_perp/Cs0 MODEL.nu = 0.01; % collisionality nu*L_perp/Cs0
% temperature ratio T_a/T_e % temperature ratio T_a/T_e
MODEL.tau_e = 1.0; MODEL.tau_e = 1.0;
MODEL.tau_i = 1.0; MODEL.tau_i = 1.0;
...@@ -36,16 +36,16 @@ MODEL.sigma_i = 1.0; ...@@ -36,16 +36,16 @@ MODEL.sigma_i = 1.0;
MODEL.q_e =-1.0; MODEL.q_e =-1.0;
MODEL.q_i = 1.0; MODEL.q_i = 1.0;
% gradients L_perp/L_x % gradients L_perp/L_x
MODEL.eta_n = 0.0; % Density MODEL.eta_n = 1.0; % Density
MODEL.eta_T = 0.0; % Temperature MODEL.eta_T = 0.0; % Temperature
MODEL.eta_B = 0.0; % Magnetic MODEL.eta_B = 0.5; % Magnetic
% Debye length % Debye length
MODEL.lambdaD = 0.0; MODEL.lambdaD = 0.0;
%% Time integration and intialization parameters %% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4'''; TIME_INTEGRATION.numerical_scheme = '''RK4''';
BASIC.nrun = 100000; BASIC.nrun = 100000;
BASIC.dt = 0.01; BASIC.dt = 0.05;
BASIC.tmax = 5.0; BASIC.tmax = 100.0;
INITIAL.initback_moments = 0.01; INITIAL.initback_moments = 0.01;
INITIAL.initnoise_moments = 0.; INITIAL.initnoise_moments = 0.;
INITIAL.iseed = 42; INITIAL.iseed = 42;
...@@ -85,92 +85,46 @@ end ...@@ -85,92 +85,46 @@ end
default_plots_options default_plots_options
SAVEFIG = 1; SAVEFIG = 1;
filename = 'results_00.h5'; filename = 'results_00.h5';
TITLE = []; default_plots_options
TITLE = [TITLE,'$\eta_n=',num2str(MODEL.eta_n),'$, ']; TITLE = [', $k_z=',num2str(GRID.kzmin)];
TITLE = [TITLE,'$\eta_B=',num2str(MODEL.eta_B),'$, '];
TITLE = [TITLE,'$\eta_T=',num2str(MODEL.eta_T),'$, '];
TITLE = [TITLE, '$\nu=',num2str(MODEL.nu),'$, '];
TITLE = [TITLE, '$(P,J)=(',num2str(GRID.pmaxe),',',num2str(GRID.jmaxe),')$'];
if MODEL.CO == -1; CONAME = 'FC';
elseif MODEL.CO == -2; CONAME = 'DC';
elseif MODEL.CO == 0; CONAME = 'LB'; end;
bare = @(pp,jj) (GRID.jmaxe +1)*pp + jj+1; % electron 1D-index
bari = @(pp,jj) bare(GRID.pmaxe, GRID.jmaxe)+(GRID.jmaxi +1)*pp + jj+1; % ion 1D-index
%% Nepj
%
% moment = 'moments_e';
%
% kr = h5read(filename,['/data/var5d/' moment '/coordkr']);
% kz = h5read(filename,['/data/var5d/' moment '/coordkz']);
% time = h5read(filename,'/data/var5d/time');
% Nepj = zeros(GRID.pmaxe+1, GRID.jmaxe+1,numel(kr),numel(kz),numel(time));
% for it = 1:numel(time)
% tmp = h5read(filename,['/data/var5d/', moment,'/', num2str(it,'%06d')]);
% Nepj(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
% end
%
% fig = figure;
%
% %HeLaZ results
% x1 = time;
% y1 = x1;
% ic = 1;
% for ikr = 1:GRID.nkr
% for ikz = 1:GRID.nkz
% for ip = 0:1
% for ij = 0:1
% for it = 1:numel(time)
% y1(it) = abs(Nepj(ip+1,ij+1,ikr,ikz,it));
% end
% semilogy(x1,y1,'-','DisplayName',...
% ['HeLaZ $N_e^{',num2str(ip),num2str(ij),'}$'],...
% 'color', line_colors(ic,:))
% hold on
% ic = ic + 1;
% end
% end
% end
% end
%
% %MOLI results
% x2 = results.time;
% y2 = x2;
% ic = 1;
% for ip = 0:1
% for ij = 0:1
% for it = 1:numel(x2)
% y2(it) = abs(results.Napj(it,bare(ip,ij)));
% end
% semilogy(x2,y2,'--','DisplayName',...
% ['MOLI $N_e^{',num2str(ip),num2str(ij),'}$'],...
% 'color', line_colors(ic,:))
% hold on
% ic = ic + 1;
% end
% end
%
% title(TITLE);
% grid on
% legend('show')
% xlabel('$t$')
% ylabel(['$|N_e^{pj}|$'])
% if SAVEFIG; FIGNAME = ['Nepj_kz_',num2str(GRID.kzmin)]; save_figure; end;
%
%% Nipj %% Nipj
% load HeLaZ data
moment = 'Ni00';
kr = h5read(filename,['/data/var2d/' moment '/coordkr']);
kz = h5read(filename,['/data/var2d/' moment '/coordkz']);
timeNi = squeeze(h5read(filename,'/data/var2d/time'));
Ni00 = zeros(numel(kr),numel(kz),numel(timeNi));
for it = 1:numel(timeNi)
tmp = h5read(filename,['/data/var2d/', moment,'/', num2str(it,'%06d')]);
Ni00(:,:,it) = tmp.real + 1i * tmp.imaginary;
end
moment = 'moments_i'; moment = 'moments_i';
kr = h5read(filename,['/data/var5d/' moment '/coordkr']); kr = h5read(filename,['/data/var5d/' moment '/coordkr']);
kz = h5read(filename,['/data/var5d/' moment '/coordkz']); kz = h5read(filename,['/data/var5d/' moment '/coordkz']);
time = h5read(filename,'/data/var5d/time'); time = h5read(filename,'/data/var5d/time');
Nipj = zeros(GRID.pmaxe+1, GRID.jmaxe+1,numel(kr),numel(kz),numel(time)); Nipj = zeros(GRID.pmaxi+1, GRID.jmaxi+1,numel(kr),numel(kz),numel(time));
for it = 1:numel(time) for it = 1:numel(time)
tmp = h5read(filename,['/data/var5d/', moment,'/', num2str(it,'%06d')]); tmp = h5read(filename,['/data/var5d/', moment,'/', num2str(it,'%06d')]);
Nipj(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary; Nipj(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
end end
err_ = mean(abs(squeeze(Nipj(1,1,1,1,:))-squeeze((Ni00(1,1,:)))));
disp(['2D 5D error :',num2str(err_)])
% growth rate :
gammaHeLaZ = LinearFit_s(squeeze(time),...
squeeze(abs(Nipj(1,1,1,1,:)))) ;
gammaMOLI = LinearFit(results.Napj,results.time,params,options);
disp(['Growth rate HeLaZ : ',num2str(gammaHeLaZ)])
disp(['Growth rate MOLI : ',num2str( gammaMOLI)])
disp(['Error : ',num2str( abs(gammaMOLI-gammaHeLaZ))])
% Plot moments
fig = figure; fig = figure;
x1 = time; x1 = time;
...@@ -284,11 +238,34 @@ plot(timephiMOLI,ERR2,'--','DisplayName','Imag') ...@@ -284,11 +238,34 @@ plot(timephiMOLI,ERR2,'--','DisplayName','Imag')
%title(TITLE); %title(TITLE);
grid on grid on
xlabel('$t$') xlabel('$t$')
ylabel('$|e_\phi|$') ylabel('$e(\phi)$')
legend('show') legend('show')
%if SAVEFIG; FIGNAME = ['err_phi_kz_',num2str(GRID.kzmin)]; save_figure; end; %if SAVEFIG; FIGNAME = ['err_phi_kz_',num2str(GRID.kzmin)]; save_figure; end;
%% Growth rate fit quantity (Ni00)
subplot(326);
y1 = squeeze(abs(Nipj(1,1,1,1,:))); % HeLaZ
y2 = squeeze(abs(results.Napj(:,bari(0,0))));
ERR3 = abs(y2-y1);
yyaxis left
semilogy(x1,y1,'-','DisplayName','HeLaZ',...
'color', line_colors(1,:)); hold on;
semilogy(x2,y2,'--','DisplayName','MOLI',...
'color', 'k')
grid on
xlabel('$t$')
ylabel('$|N_i^{00}|$')
legend('show')
yyaxis right
semilogy(x1, ERR3,'color', line_colors(2,:));
grid on
xlabel('$t$')
ylabel('$e(N_i^{00})$')
%% Finish and save
suptitle(TITLE); suptitle(TITLE);
if SAVEFIG; FIGNAME = ['kz_',num2str(GRID.kzmin)]; save_figure; end; if SAVEFIG; FIGNAME = ['kz_',num2str(GRID.kzmin)]; save_figure; end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
wk/fort.90
+ 16
− 16
View file @ ef30a980
&BASIC &BASIC
nrun=100000 nrun=100000
dt=0.01 dt=0.05
tmax=5 ! time normalized to 1/omega_pe tmax=100 ! time normalized to 1/omega_pe
/ /
&GRID &GRID
pmaxe =6 ! Electron Hermite moments pmaxe =12 ! Electron Hermite moments
jmaxe = 6 ! Electron Laguerre moments jmaxe = 7 ! Electron Laguerre moments
pmaxi = 6 ! Ion Hermite moments pmaxi = 12 ! Ion Hermite moments
jmaxi = 6 ! Ion Laguerre moments jmaxi = 7 ! Ion Laguerre moments
nkr = 1 nkr = 1
krmin = 0 krmin = 0
krmax = 0 krmax = 0
nkz = 1 nkz = 20
kzmin = 1 kzmin = 0.1
kzmax = 1 kzmax = 2
/ /
&OUTPUT_PAR &OUTPUT_PAR
nsave_0d = 0 nsave_0d = 0
nsave_1d = 0 nsave_1d = 0
nsave_2d = 1 nsave_2d = 1
nsave_5d = 1 nsave_5d = 0
write_Ni00 = .true. write_Ni00 = .true.
write_moments = .true. write_moments = .true.
write_phi = .true. write_phi = .true.
...@@ -29,16 +29,16 @@ ...@@ -29,16 +29,16 @@
&MODEL_PAR &MODEL_PAR
! Collisionality ! Collisionality
CO = -2 ! Collision operator (-1:Full Coulomb, 0: Dougherty) CO = -2 ! Collision operator (-1:Full Coulomb, 0: Dougherty)
nu = 1 ! Normalized collision frequency normalized to plasma frequency nu = 0.01 ! Normalized collision frequency normalized to plasma frequency
tau_e = 1 ! T_e/T_e tau_e = 1 ! T_e/T_e
tau_i = 1 ! T_i/T_e temperature ratio tau_i = 1 ! T_i/T_e temperature ratio
sigma_e = 0.023338 ! sqrt(m_e/m_i) mass ratio sigma_e = 0.023338 ! sqrt(m_e/m_i) mass ratio
sigma_i = 1 ! sqrt(m_i/m_i) sigma_i = 1 ! sqrt(m_i/m_i)
q_e = -1 ! Electrons charge q_e = -1 ! Electrons charge
q_i = 1 ! Ions charge q_i = 1 ! Ions charge
eta_n = 0 ! L_perp / L_n Density gradient eta_n = 1 ! L_perp / L_n Density gradient
eta_T = 0 ! L_perp / L_T Temperature gradient eta_T = 0 ! L_perp / L_T Temperature gradient
eta_B = 0 ! L_perp / L_B Magnetic gradient and curvature eta_B = 0.5 ! L_perp / L_B Magnetic gradient and curvature
lambdaD = 0 ! Debye length lambdaD = 0 ! Debye length
/ /
&INITIAL_CON &INITIAL_CON
...@@ -47,9 +47,9 @@ ...@@ -47,9 +47,9 @@
! Noise amplitude ! Noise amplitude
initnoise_moments=0 initnoise_moments=0
iseed=42 iseed=42
selfmat_file='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_6_Jmaxe_6_Pmaxi_6_Jmaxi_6_pamaxx_10.h5' selfmat_file='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_12_Jmaxe_7_Pmaxi_12_Jmaxi_7_pamaxx_10.h5'
eimat_file='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_6_Jmaxe_6_Pmaxi_6_Jmaxi_6_pamaxx_10_tau_1.0000_mu_0.0233.h5' eimat_file='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_12_Jmaxe_7_Pmaxi_12_Jmaxi_7_pamaxx_10_tau_1.0000_mu_0.0233.h5'
iemat_file='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_6_Jmaxe_6_Pmaxi_6_Jmaxi_6_pamaxx_10_tau_1.0000_mu_0.0233.h5' iemat_file='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_12_Jmaxe_7_Pmaxi_12_Jmaxi_7_pamaxx_10_tau_1.0000_mu_0.0233.h5'
/ /
&TIME_INTEGRATION_PAR &TIME_INTEGRATION_PAR
numerical_scheme='RK4' numerical_scheme='RK4'
......
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