From ef30a9806ed826a5a7ef5e04b8d892e0ca656c5a Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 7 Jul 2020 10:12:20 +0200
Subject: [PATCH] Cleaned and refreshed
---
wk/benchmark_kperp_scan.m | 80 +++++++++++----------
wk/benchmark_time_solver.m | 141 ++++++++++++++++---------------------
wk/fort.90 | 32 ++++-----
3 files changed, 118 insertions(+), 135 deletions(-)
diff --git a/wk/benchmark_kperp_scan.m b/wk/benchmark_kperp_scan.m
index bdec0eb9..e14e54eb 100644
--- a/wk/benchmark_kperp_scan.m
+++ b/wk/benchmark_kperp_scan.m
@@ -12,19 +12,19 @@ OUTPUTS.write_phi = '.true.';
OUTPUTS.write_doubleprecision = '.true.';
OUTPUTS.resfile0 = '''results''';
%% Grid parameters
-GRID.pmaxe = 15;
-GRID.jmaxe = 6;
-GRID.pmaxi = 15;
-GRID.jmaxi = 6;
+GRID.pmaxe = 12;
+GRID.jmaxe = 7;
+GRID.pmaxi = 12;
+GRID.jmaxi = 7;
GRID.nkr = 1;
GRID.krmin = 0.;
GRID.krmax = 0.;
GRID.nkz = 20;
GRID.kzmin = 0.1;
-GRID.kzmax = 1.5;
+GRID.kzmax = 2.0;
%% Model parameters
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
MODEL.tau_e = 1.0;
MODEL.tau_i = 1.0;
@@ -48,7 +48,18 @@ BASIC.tmax = 100.0;
INITIAL.initback_moments = 0.01;
INITIAL.initnoise_moments = 0.;
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
INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
@@ -72,63 +83,58 @@ end
%% Analysis and basic figures
SAVEFIG = 1;
filename = 'results_00.h5';
-
-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),')$'];
+default_plots_options
%% Growth rate analysis
-gammas = zeros(numel(kr),numel(kz));
-shifts = zeros(numel(kr),numel(kz));
-
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));
+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')]);
- Nipj(it,:,:) = tmp.real + 1i * tmp.imaginary;
+ Ni00(:,:,it) = tmp.real + 1i * tmp.imaginary;
end
-% Linear fit of log(Napj)
-x1 = timeNi;
-itmin = ceil(0.5 * numel(timeNi)); %Take the second half of the time evolution
+% Amplitude ratio method
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)
- fit = polyfit(x1(itmin:end),log(abs(Nipj(itmin:end,ikr,ikz))),1);
- gammas(ikr,ikz) = fit(1);
- shifts(ikr,ikz) = fit(2);
+ gammas(ikr,ikz) = LinearFit_s(squeeze(timeNi), squeeze(abs(Ni00(ikr,ikz,:))));
end
%% Plot
fig = figure;
+X = kz*sqrt(1+MODEL.tau_i); % Convert to Ricci 2006 Normalization
+
+subplot(121) % growth rate
%HeLaZ results
-X = kz*sqrt(1+MODEL.tau_i);
-Y = gammas(1,:);
-plot(X,Y,'-','DisplayName','HeLaZ')
+Y1 = gammas(1,:);
+plot(X,Y1,'-','DisplayName','HeLaZ')
hold on
%MOLI results
-X = kz*sqrt(1+MODEL.tau_i);
-Y = results.kperp.Maxgammas;
-plot(X,Y,'--','DisplayName','MOLI');
+Y2 = results.kperp.Maxgammas;
+plot(X,Y2,'--','DisplayName','MOLI');
title(TITLE);
grid on
legend('show')
xlabel('$k_\perp * (1+\tau)^{1/2}$')
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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
diff --git a/wk/benchmark_time_solver.m b/wk/benchmark_time_solver.m
index a330b5f3..2eaaaa26 100644
--- a/wk/benchmark_time_solver.m
+++ b/wk/benchmark_time_solver.m
@@ -21,11 +21,11 @@ GRID.nkr = 1;
GRID.krmin = 0.;
GRID.krmax = 0.;
GRID.nkz = 1;
-GRID.kzmin = 1.0;
-GRID.kzmax = 1.0;
+GRID.kzmin = 0.67;
+GRID.kzmax = 0.67;
%% Model parameters
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
MODEL.tau_e = 1.0;
MODEL.tau_i = 1.0;
@@ -36,16 +36,16 @@ MODEL.sigma_i = 1.0;
MODEL.q_e =-1.0;
MODEL.q_i = 1.0;
% 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_B = 0.0; % Magnetic
+MODEL.eta_B = 0.5; % Magnetic
% Debye length
MODEL.lambdaD = 0.0;
%% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4''';
BASIC.nrun = 100000;
-BASIC.dt = 0.01;
-BASIC.tmax = 5.0;
+BASIC.dt = 0.05;
+BASIC.tmax = 100.0;
INITIAL.initback_moments = 0.01;
INITIAL.initnoise_moments = 0.;
INITIAL.iseed = 42;
@@ -85,92 +85,46 @@ end
default_plots_options
SAVEFIG = 1;
filename = 'results_00.h5';
-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),')$'];
-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
+default_plots_options
+TITLE = [', $k_z=',num2str(GRID.kzmin)];
-%% 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
+% 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';
kr = h5read(filename,['/data/var5d/' moment '/coordkr']);
kz = h5read(filename,['/data/var5d/' moment '/coordkz']);
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)
tmp = h5read(filename,['/data/var5d/', moment,'/', num2str(it,'%06d')]);
Nipj(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
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;
x1 = time;
@@ -284,11 +238,34 @@ plot(timephiMOLI,ERR2,'--','DisplayName','Imag')
%title(TITLE);
grid on
xlabel('$t$')
-ylabel('$|e_\phi|$')
+ylabel('$e(\phi)$')
legend('show')
%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);
if SAVEFIG; FIGNAME = ['kz_',num2str(GRID.kzmin)]; save_figure; end;
-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
diff --git a/wk/fort.90 b/wk/fort.90
index 3c697bf2..2017fc90 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -1,25 +1,25 @@
&BASIC
nrun=100000
- dt=0.01
- tmax=5 ! time normalized to 1/omega_pe
+ dt=0.05
+ tmax=100 ! time normalized to 1/omega_pe
/
&GRID
- pmaxe =6 ! Electron Hermite moments
- jmaxe = 6 ! Electron Laguerre moments
- pmaxi = 6 ! Ion Hermite moments
- jmaxi = 6 ! Ion Laguerre moments
+ pmaxe =12 ! Electron Hermite moments
+ jmaxe = 7 ! Electron Laguerre moments
+ pmaxi = 12 ! Ion Hermite moments
+ jmaxi = 7 ! Ion Laguerre moments
nkr = 1
krmin = 0
krmax = 0
- nkz = 1
- kzmin = 1
- kzmax = 1
+ nkz = 20
+ kzmin = 0.1
+ kzmax = 2
/
&OUTPUT_PAR
nsave_0d = 0
nsave_1d = 0
nsave_2d = 1
- nsave_5d = 1
+ nsave_5d = 0
write_Ni00 = .true.
write_moments = .true.
write_phi = .true.
@@ -29,16 +29,16 @@
&MODEL_PAR
! Collisionality
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_i = 1 ! T_i/T_e temperature ratio
sigma_e = 0.023338 ! sqrt(m_e/m_i) mass ratio
sigma_i = 1 ! sqrt(m_i/m_i)
q_e = -1 ! Electrons 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_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
/
&INITIAL_CON
@@ -47,9 +47,9 @@
! Noise amplitude
initnoise_moments=0
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'
- 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'
- 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'
+ 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_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_12_Jmaxe_7_Pmaxi_12_Jmaxi_7_pamaxx_10_tau_1.0000_mu_0.0233.h5'
/
&TIME_INTEGRATION_PAR
numerical_scheme='RK4'
--
GitLab