From bb7b5ed9b9f19bfddf364233621eea2637e4cb97 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 15 Jun 2021 14:06:54 +0200
Subject: [PATCH] scripts updates

---
 matlab/load_results.m              |  12 ++
 matlab/photomaton.m                |  14 +-
 matlab/setup.m                     |   6 +-
 matlab/write_fort90.m              |   2 +
 wk/analysis_2D.m                   |  90 ++++++++++---
 wk/linear_study.m                  |   4 +-
 wk/load_multiple_outputs_marconi.m |   3 +-
 wk/local_run.m                     |  20 +--
 wk/marconi_run.m                   |  16 +--
 wk/open_figure_script.m            |  17 ++-
 wk/transport_results_2_6.m         | 208 ++++++++++++++---------------
 11 files changed, 227 insertions(+), 165 deletions(-)

diff --git a/matlab/load_results.m b/matlab/load_results.m
index c4489bc0..412b92a4 100644
--- a/matlab/load_results.m
+++ b/matlab/load_results.m
@@ -11,6 +11,8 @@ W_PHI     = strcmp(h5readatt(filename,'/data/input','write_phi')  ,'y');
 W_NA00    = strcmp(h5readatt(filename,'/data/input','write_Na00') ,'y');
 W_NAPJ    = strcmp(h5readatt(filename,'/data/input','write_Napj') ,'y');
 W_SAPJ    = strcmp(h5readatt(filename,'/data/input','write_Sapj') ,'y');
+W_DENS    = strcmp(h5readatt(filename,'/data/input','write_dens') ,'y');
+W_TEMP    = strcmp(h5readatt(filename,'/data/input','write_temp') ,'y');
 
 
 if W_GAMMA
@@ -37,3 +39,13 @@ if W_SAPJ
     [Sipj, Ts5D, dt5D] = load_5D_data(filename, 'Sipj');
      Sepj              = load_5D_data(filename, 'Sepj');
 end
+
+if W_DENS
+    [DENS_E, Ts2D, dt2D] = load_2D_data(filename, 'dens_e');
+    [DENS_I, Ts2D, dt2D] = load_2D_data(filename, 'dens_i');
+end
+
+if W_TEMP
+    [TEMP_E, Ts2D, dt2D] = load_2D_data(filename, 'temp_e');
+    [TEMP_I, Ts2D, dt2D] = load_2D_data(filename, 'temp_i');
+end
\ No newline at end of file
diff --git a/matlab/photomaton.m b/matlab/photomaton.m
index d5508bfc..1d3d0db4 100644
--- a/matlab/photomaton.m
+++ b/matlab/photomaton.m
@@ -1,13 +1,13 @@
 if 0
 %% Photomaton : real space
-% FIELD = ni00; FNAME = 'ni'; FIELDLTX = '$N_i^{00}$'; XX = RR; YY = ZZ;
+FIELD = ni00; FNAME = 'ni'; FIELDLTX = '$N_i^{00}$'; XX = RR; YY = ZZ;
 % FIELD = ne00; FNAME = 'ne'; FIELDLTX = '$N_e^{00}$'; XX = RR; YY = ZZ;
-FIELD = phi; FNAME = 'phi'; FIELDLTX = '$\phi$'; XX = RR; YY = ZZ;
-% FIELD = dr2phi; FNAME = 'dr2phi'; XX = RR; YY = ZZ;
-tf = 100;  [~,it1] = min(abs(Ts2D-tf));
-tf = 285;  [~,it2] = min(abs(Ts2D-tf)); 
-tf = 350; [~,it3] = min(abs(Ts2D-tf));
-tf = 1100; [~,it4] = min(abs(Ts2D-tf));
+% FIELD = phi; FNAME = 'phi'; FIELDLTX = '$\phi$'; XX = RR; YY = ZZ;
+% FIELD = dr2phi; FNAME = 'dr2phi'; FIELDLTX = '$\partial_r^2\phi$'; XX = RR; YY = ZZ;
+tf = 10;  [~,it1] = min(abs(Ts2D-tf));
+tf = 40;  [~,it2] = min(abs(Ts2D-tf)); 
+tf = 80; [~,it3] = min(abs(Ts2D-tf));
+tf = 100; [~,it4] = min(abs(Ts2D-tf));
 fig = figure; FIGNAME = [FNAME,'_snaps','_',PARAMS]; set(gcf, 'Position',  [100, 100, 1500, 350])
 plt = @(x) x;%./max(max(x));
     subplot(141)
diff --git a/matlab/setup.m b/matlab/setup.m
index dd176526..270c8db0 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -163,6 +163,8 @@ if W_PHI;    OUTPUTS.write_phi   = '.true.'; else; OUTPUTS.write_phi   = '.false
 if W_NA00;   OUTPUTS.write_Na00  = '.true.'; else; OUTPUTS.write_Na00  = '.false.';end;
 if W_NAPJ;   OUTPUTS.write_Napj  = '.true.'; else; OUTPUTS.write_Napj  = '.false.';end;
 if W_SAPJ;   OUTPUTS.write_Sapj  = '.true.'; else; OUTPUTS.write_Sapj  = '.false.';end;
+if W_DENS;   OUTPUTS.write_dens  = '.true.'; else; OUTPUTS.write_dens  = '.false.';end;
+if W_TEMP;   OUTPUTS.write_temp  = '.true.'; else; OUTPUTS.write_temp  = '.false.';end;
 OUTPUTS.resfile0    = '''outputs''';
 OUTPUTS.rstfile0    = '''checkpoint''';
 OUTPUTS.job2load    = JOB2LOAD;
@@ -176,8 +178,8 @@ end
 %% Compile and WRITE input file
 INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
 nproc = 1;
-% MAKE  = 'cd ..; make; cd wk';
-% system(MAKE);
+MAKE  = 'cd ..; make; cd wk';
+system(MAKE);
 %%
 disp(['Set up ',SIMID]);
 disp([resolution,gridname,degngrad]);
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index ed9b721d..97dc7b07 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -35,6 +35,8 @@ fprintf(fid,['  write_phi   = ', OUTPUTS.write_phi,'\n']);
 fprintf(fid,['  write_Na00 = ', OUTPUTS.write_Na00,'\n']);
 fprintf(fid,['  write_Napj = ', OUTPUTS.write_Napj,'\n']);
 fprintf(fid,['  write_Sapj = ', OUTPUTS.write_Sapj,'\n']);
+fprintf(fid,['  write_dens = ', OUTPUTS.write_dens,'\n']);
+fprintf(fid,['  write_temp = ', OUTPUTS.write_temp,'\n']);
 fprintf(fid,['  resfile0      = ', OUTPUTS.resfile0,'\n']);
 fprintf(fid,['  rstfile0      = ', OUTPUTS.rstfile0,'\n']);
 fprintf(fid,['  job2load      = ', num2str(OUTPUTS.job2load),'\n']);
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index ee087aa6..6eebb675 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -5,18 +5,17 @@ if 1 % Local results
 outfile ='';
 outfile ='';
 outfile ='';
-outfile ='kobayashi/100x50_L_50_P_4_J_2_eta_0.71429_nu_5e-03_DGGK_CLOS_0_mu_1e-02';
-% outfile ='kobayashi/100x50_L_50_P_4_J_2_eta_0.71429_nu_5e-03_SGGK_CLOS_0_mu_1e-02';
-% outfile ='v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.5_nu_1e-01_DGGK_CLOS_0_mu_2e-02';
-
+outfile ='test_diagnostics/100x50_L_50_P_2_J_1_eta_0.5_nu_5e-03_DGGK_CLOS_0_mu_1e-02';
     BASIC.RESDIR      = ['../results/',outfile,'/'];
     BASIC.MISCDIR     = ['../results/',outfile,'/'];
 end
 if 0 % Marconi results
 outfile ='';
 outfile ='';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/kobayashi/100x50_L_50_P_6_J_3_eta_0.71429_nu_5e-03_DGGK_CLOS_0_mu_1e-02/out.txt';
-% outfile =['/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_',num2str(ETA_),'_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt'];
+outfile ='';
+outfile ='';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt';
 % load_marconi(outfile);
 BASIC.RESDIR      = ['../',outfile(46:end-8),'/'];
 BASIC.MISCDIR     = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
@@ -24,7 +23,7 @@ end
 
 %%
 % JOBNUM = 1; load_results;
-JOBNUMMAX = 0; compile_results %Compile the results from first output found to JOBNUMMAX if existing
+JOBNUMMAX = 20; compile_results %Compile the results from first output found to JOBNUMMAX if existing
 
 %% Retrieving max polynomial degree and sampling info
 Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe);
@@ -60,16 +59,26 @@ dzni00 = zeros(Nr,Nz,Ns2D);
 np_i   = zeros(Nr,Nz,Ns5D); % Ion particle density
 si00   = zeros(Nr,Nz,Ns5D);
 phi    = zeros(Nr,Nz,Ns2D);
+dens_e = zeros(Nr,Nz,Ns2D);
+dens_i = zeros(Nr,Nz,Ns2D);
+temp_e = zeros(Nr,Nz,Ns2D);
+temp_i = zeros(Nr,Nz,Ns2D);
 drphi  = zeros(Nr,Nz,Ns2D);
 dzphi  = zeros(Nr,Nz,Ns2D);
 dr2phi = zeros(Nr,Nz,Ns2D);
 
 for it = 1:numel(Ts2D)
     NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it);
-    ne00(:,:,it)  = real(fftshift(ifft2((NE_),Nr,Nz)));
-    ni00(:,:,it)  = real(fftshift(ifft2((NI_),Nr,Nz)));
-    dzni00(:,:,it) = real(fftshift(ifft2(1i*KZ.*(NI_),Nr,Nz)));
-    phi (:,:,it)  = real(fftshift(ifft2((PH_),Nr,Nz)));
+    DENS_E_ = DENS_E(:,:,it); DENS_I_ = DENS_I(:,:,it);
+    TEMP_E_ = TEMP_E(:,:,it); TEMP_I_ = TEMP_I(:,:,it);
+    ne00(:,:,it)    = real(fftshift(ifft2((NE_),Nr,Nz)));
+    ni00(:,:,it)    = real(fftshift(ifft2((NI_),Nr,Nz)));
+    dzni00(:,:,it)  = real(fftshift(ifft2(1i*KZ.*(NI_),Nr,Nz)));
+    phi (:,:,it)    = real(fftshift(ifft2((PH_),Nr,Nz)));
+    dens_e (:,:,it) = real(fftshift(ifft2((DENS_E_),Nr,Nz)));
+    dens_i (:,:,it) = real(fftshift(ifft2((DENS_I_),Nr,Nz)));
+    temp_e (:,:,it) = real(fftshift(ifft2((TEMP_E_),Nr,Nz)));
+    temp_i (:,:,it) = real(fftshift(ifft2((TEMP_I_),Nr,Nz)));
     drphi(:,:,it) = real(fftshift(ifft2(1i*KR.*(PH_),Nr,Nz)));
     dr2phi(:,:,it)= real(fftshift(ifft2(-KR.^2.*(PH_),Nr,Nz)));
     dzphi(:,:,it) = real(fftshift(ifft2(1i*KZ.*(PH_),Nr,Nz)));
@@ -329,7 +338,7 @@ fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position',  [100, 10
 save_figure
 end
 
-if 1
+if 0
 %% Averaged shear and Reynold stress profiles
 figure;
 plt = @(x) squeeze(mean(mean(real(x(:,:,its2D:ite2D)),2),3))/max(abs(squeeze(mean(mean(real(x(:,:,its2D:ite2D)),2),3))));
@@ -340,8 +349,37 @@ plot(r,plt(-drphi.*dzphi-drphi.*dzni00),'--','DisplayName','$\Pi_\phi+\Pi_{ni00}
 xlim([-60,60]); xlabel('$r/\rho_s$');
 end
 
+if 0
+%% |phi_k|^2 spectra (Kobayashi 2015 fig 3)
+tstart = 500; tend = 5000;
+[~,itstart] = min(abs(Ts2D-tstart));
+[~,itend]   = min(abs(Ts2D-tend));
+trange = itstart:itend;
+%full kperp points
+phi_k_2 = reshape(mean((abs(PHI(:,:,trange))).^2,3),[numel(KR),1]);
+kperp  = reshape(sqrt(KR.^2+KZ.^2),[numel(KR),1]);
+%interpolated kperps
+nk_noAA = floor(2/3*numel(kr));
+kp_ip = kr(1:nk_noAA);
+[thg, rg] = meshgrid(linspace(-pi/2,pi/2,2*nk_noAA),kp_ip);
+[xn,yn] = pol2cart(thg,rg);
+[xc,yc] = meshgrid(sort(kz),kr);
+Z_rth = interp2(xc,yc,mean((abs(PHI(:,:,trange))).^2,3),xn,yn);
+z_k = mean(Z_rth,2);
+%for theorical trends
+kp_1d  = linspace(0.1,5,100);
+figure;
+scatter(kperp,phi_k_2,'.'); hold on; grid on;
+plot(kp_ip,z_k);
+plot(kp_1d,1e4*kp_1d.^(-13/3));
+plot(kp_1d,1e4*kp_1d.^(-3)./(1+kp_1d.^2).^2);
+set(gca, 'XScale', 'log');set(gca, 'YScale', 'log');
+xlabel('$k_\perp \rho_s$'); ylabel('$|\phi_k|^2$');
+xlim([0.1,10]);
+end
+
 %%
-t0    = 00;
+t0    = 000;
 [~, it02D] = min(abs(Ts2D-t0));
 [~, it05D] = min(abs(Ts5D-t0));
 skip_ = 1; 
@@ -350,18 +388,34 @@ FRAMES_2D = it02D:skip_:numel(Ts2D);
 FRAMES_5D = it05D:skip_:numel(Ts5D);
 %% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 if 0
-%% Density ion
+%% part density ion
 GIFNAME = ['ni',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
-FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
+FIELD = real(dens_i); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
 FIELDNAME = '$n_i$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
 % create_gif
 create_mov
 end
 if 0
-%% Density electrons
-GIFNAME = ['ne',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
+%% part temperature ion
+GIFNAME = ['Ti',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
+FIELD = real(temp_i); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
+FIELDNAME = '$T_i$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
+% create_gif
+create_mov
+end
+if 0
+%% GC Density ion
+GIFNAME = ['ni00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
+FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
+FIELDNAME = '$n_i^{00}$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
+% create_gif
+create_mov
+end
+if 0
+%% GC Density electrons
+GIFNAME = ['ne00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
 FIELD = real(ne00); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
-FIELDNAME = '$n_e$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
+FIELDNAME = '$n_e^{00}$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
 % create_gif
 create_mov
 end
diff --git a/wk/linear_study.m b/wk/linear_study.m
index 671fa975..58855861 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -18,8 +18,8 @@ NU_HYP  = 0.0;   % Hyperdiffusivity coefficient
 LAMBDAD = 0.0;
 NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
-N       = 66;     % Frequency gridpoints (Nkr = N/2)
-L       = 150;     % Size of the squared frequency domain
+N       = 100;     % Frequency gridpoints (Nkr = N/2)
+L       = 100;     % Size of the squared frequency domain
 KREQ0   = 1;      % put kr = 0
 MU_P    = 0.0;     % Hermite  hyperdiffusivity -mu_p*(d/dvpar)^4 f
 MU_J    = 0.0;     % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
diff --git a/wk/load_multiple_outputs_marconi.m b/wk/load_multiple_outputs_marconi.m
index 16b333b5..0c7c81ea 100644
--- a/wk/load_multiple_outputs_marconi.m
+++ b/wk/load_multiple_outputs_marconi.m
@@ -1,7 +1,8 @@
 addpath(genpath('../matlab')) % ... add
 %% Paste the list of simulation results to load
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.5_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
 load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
 load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.7_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
 load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.8_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
 load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.9_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.5_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/kobayashi/100x50_L_50_P_10_J_5_eta_0.71429_nu_5e-03_DGGK_CLOS_0_mu_1e-02/out.txt')
diff --git a/wk/local_run.m b/wk/local_run.m
index 8d3c7ac3..7d099685 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -5,19 +5,19 @@ CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
 NU      = 5.0e-3;   % Collision frequency
-ETAB    = 1/1.4;    % Magnetic gradient
+ETAB    = 0.5;%1/1.4;    % Magnetic gradient
 ETAN    = 1.0;    % Density gradient
 NU_HYP  = 1.0;
 %% GRID PARAMETERS
 N       = 100;     % Frequency gridpoints (Nkr = N/2)
 L       = 50;     % Size of the squared frequency domain
-P       = 4;
-J       = 2;
+P       = 2;
+J       = 1;
 MU_P    = 0.0;     % Hermite  hyperdiffusivity -mu_p*(d/dvpar)^4 f
 MU_J    = 0.0;     % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
 %% TIME PARAMETERS
-TMAX    = 2000;  % Maximal time unit
-DT      = 1e-2;   % Time step
+TMAX    = 500;  % Maximal time unit
+DT      = 2e-2;   % Time step
 SPS0D   = 1;    % Sampling per time unit for profiler
 SPS2D   = 1;      % Sampling per time unit for 2D arrays
 SPS5D   = 1/50;    % Sampling per time unit for 5D arrays
@@ -27,12 +27,12 @@ JOB2LOAD= 0;
 %% OPTIONS AND NAMING
 % Collision operator
 % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
-CO      = 2;
+CO      =1;
 CLOS    = 0;   % Closure model (0: =0 truncation)
 NL_CLOS = -1;   % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
-% SIMID   = 'test';  % Name of the simulation
-SIMID   = 'kobayashi';  % Name of the simulation
-% SIMID   = ['v2.6_P_',num2str(P),'_J_',num2str(J)];  % Name of the simulation
+SIMID   = 'test_diagnostics';  % Name of the simulation
+% SIMID   = 'kobayashi';  % Name of the simulation
+% SIMID   = ['v2.7_P_',num2str(P),'_J_',num2str(J)];  % Name of the simulation
 NON_LIN = 1;   % activate non-linearity (is cancelled if KREQ0 = 1)
 INIT_ZF = 0; ZF_AMP = 0.0;
 %% OUTPUTS
@@ -42,6 +42,8 @@ W_PHI    = 1;
 W_NA00   = 1;
 W_NAPJ   = 1;
 W_SAPJ   = 0;
+W_DENS   = 1;
+W_TEMP   = 1;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% unused
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 92dda8fe..eab3af94 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -1,9 +1,9 @@
 clear all;
 addpath(genpath('../matlab')) % ... add
-SUBMIT = 0; % To submit the job automatically
+SUBMIT = 1; % To submit the job automatically
 % EXECNAME = 'helaz_dbg';
-  EXECNAME = 'helaz_2.62';
-for ETAB = [0.5]
+  EXECNAME = 'helaz_2.71';
+for ETAB = [0.6]
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Set Up parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -28,22 +28,22 @@ INIT_ZF = 0; ZF_AMP = 0.0;
 N       = 200;    % Frequency gridpoints (Nkr = N/2)
 L       = 120;    % Size of the squared frequency domain
 P       = 10;     % Electron and Ion highest Hermite polynomial degree
-J       = 05;     % Electron and Ion highest Laguerre polynomial degree
+J       = 5;     % Electron and Ion highest Laguerre polynomial degree
 MU_P    = 0.0;% Hermite  hyperdiffusivity -mu_p*(d/dvpar)^4 f
 MU_J    = 0.0;% Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
 %% TIME PARAMETERS
-TMAX    = 10000;  % Maximal time unit
-DT      = 1e-3;  % Time step
+TMAX    = 5000;  % Maximal time unit
+DT      = 1e-5;  % Time step
 SPS0D   = 1;     % Sampling per time unit for profiler
 SPS2D   = 1/4;     % Sampling per time unit for 2D arrays
 SPS5D   = 1/100;  % Sampling per time unit for 5D arrays
 SPSCP   = 0;     % Sampling per time unit for checkpoints
-RESTART = 1;     % To restart from last checkpoint
+RESTART = 0;     % To restart from last checkpoint
 JOB2LOAD= 0;
 %% Naming
 % SIMID   = 'kobayashi';  % Name of the simulation
 % SIMID   = 'test';  % Name of the simulation
-SIMID   = ['v2.6_P_',num2str(P),'_J_',num2str(J)];  % Name of the simulation
+SIMID   = ['v2.7_P_',num2str(P),'_J_',num2str(J)];  % Name of the simulation
 PREFIX  =[];
 % PREFIX  = sprintf('%d_%d_',NP_P, NP_KR);
 %% Options
diff --git a/wk/open_figure_script.m b/wk/open_figure_script.m
index c69ae0c0..947dcf31 100644
--- a/wk/open_figure_script.m
+++ b/wk/open_figure_script.m
@@ -2,16 +2,15 @@ simname_ = '';
 fname_='';
 %% Marconi output file
 % fname_='';
-% fname_='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
-fname_='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.7_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt';
+% fname_='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt';
+fname_='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt';
 simname_ = fname_(54:end-8);
 
 %%
-simname_ = '';
-simname_ = '';
-simname_ = '';
-simname_ = '';
-simname_ = 'kobayashi/100x50_L_50_P_2_J_1_eta_0.71429_nu_5e-03_DGGK_CLOS_0_mu_1e-02';
+% simname_ = '';
+% simname_ = '';
+% simname_ = '';
+% simname_ = 'v2.7_P_4_J_2/100x50_L_50_P_4_J_2_eta_0.7_nu_5e-02_DGGK_CLOS_0_mu_1e-02';
 
 % simname_ = 'v2.6_P_2_J_1/200x100_L_120_P_2_J_1_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-01';
 % simname_ = 'v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02';
@@ -21,9 +20,9 @@ simname_ = 'kobayashi/100x50_L_50_P_2_J_1_eta_0.71429_nu_5e-03_DGGK_CLOS_0_mu_1e
 
 
 %%
-figname_ = '/fig/ZF_transport_drphi_';
+% figname_ = '/fig/ZF_transport_drphi_';
 % figname_ = '/fig/space_time_';
-% figname_ = '/fig/phi_shear_phase_space_';
+figname_ = '/fig/phi_shear_phase_space_';
 
 path_    = '../results/';
 
diff --git a/wk/transport_results_2_6.m b/wk/transport_results_2_6.m
index 47e78c83..b329b174 100644
--- a/wk/transport_results_2_6.m
+++ b/wk/transport_results_2_6.m
@@ -1,156 +1,146 @@
 default_plots_options
 %% nuDGGK = 0.1
 figure; set(gcf, 'Position',  [100, 100, 1200, 350])
-% 6,3 nuDGGK=1.0
-eta   = [0.5,0.6,0.7,0.8];
-gamma = [32.6443 3.6895 0.3744 0.056];
-std_g = [6.1529 0.7986 0.0493 0.0134];
+% 6,3 nuDGGK=0.1
+eta   = [0.5,0.6,0.7];
+gamma = [116 24.396 3.684];
+std_g = [23  3  0.4];
 subplot(121)
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=1.0$',...
-    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(1,:)); hold on;
+errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.1$',...
+    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(6,:)); hold on;
 subplot(122)
 
-% 10,5 nuDGGK = 1.0
-eta = [0.5 0.6,0.7,0.8];
-gamma = [32.6 3.5546 0.3917 0.0500];
-std_g = [7.7 0.5846 0.0486 0.0088];
-shear = [1.8505 0.60866 0.048249];
-std_s = [0.1599 0.00614 0.0061403];
+% 10,5 nuDGGK = 0.1
+eta = [0.5 0.6,0.7];
+gamma = [110 24.9 4.1];
+std_g = [23 4 0.38];
+shear = [4.6 4.0 1.35];
+std_s = [0.5 0.41 0.12];
 subplot(121)
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=1.0$',...
-    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(2,:)); hold on;
+errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
+    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(6,:)); hold on;
 subplot(122)
-errorbar(eta(2:end),shear,std_s,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=1.0$',...
-    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(2,:)); hold on;
+errorbar(eta,shear,std_s,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
+    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(6,:)); hold on;
 
-% 12,6 nuDGGK = 1.0
-eta = [0.6];
-gamma = [4.064];
-std_g = [0.7964];
+kmax = [0.576 0.5236 0.314];
+gmax = [0.305 0.20 0.06];
 subplot(121)
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=12,06$, $\nu_{DGGK}=1.0$',...
-    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(3,:)); hold on;
-subplot(122)
-
-
 % Mix len Ricci 2006
-subplot(121)
 semilogy([0.5  1.0],[10  1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
+% Bohm transport :
+btransp = gmax./kmax.^2;
+semilogy(.5:.1:.7,btransp','--','color',[0.4,0,0]+0.6, 'DisplayName','$\gamma/k^2$');
 set(gca, 'YScale', 'log'); grid on; legend('show')
 xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$') 
 subplot(122)
-grate = [0.2131 0.106 0.02021];
-semilogy([0.6 0.7 0.8],grate,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
+semilogy([0.5 0.6 0.7],gmax,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
 set(gca, 'YScale', 'log'); grid on; legend('show')
 xlabel('$L_n/R$'); ylabel('$V_E$''') 
 
-%% nuDGGK = 0.5
+%% nuDGGK = 0.01
 figure; set(gcf, 'Position',  [100, 100, 1200, 350])
-% 6,3 nuDGGK=0.5
+% 6,3 nuDGGK=0.01
 eta = [0.5,0.6,0.7,0.8];
-gamma = [20.511 2.6292 0.2353 0.057];
-std_g = [3.67 1.2 0.055 0.004];
-shear = [nan 1.7417  0.57345 0.25155];
-std_s = [nan 0.35991 0.041 0.00913];
+gamma = [30.2 4.3 0.37 0.06];
+std_g = [4 0.7 0.05 0.008];
+shear = [5.0 2.0 0.6 0.16];
+std_s = [0.37 0.17 0.03 0.012];
 subplot(121)
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.5$',...
-    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(4,:)); hold on;
+errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.01$',...
+    'LineWidth',2.0,'MarkerSize',5, 'Color', 'b'); hold on;
 subplot(122)
-errorbar(eta,shear,std_s,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.5$',...
-    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(4,:)); hold on;
+errorbar(eta,shear,std_s,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.01$',...
+    'LineWidth',2.0,'MarkerSize',5, 'Color', 'b'); hold on;
 
-% 10,5 nuDGGK=0.5
-eta    = [0.6,0.7,0.8 0.9];
-gamma  = [2.4949 0.23368 0.057792 0.023572];
-std_g  = [0.8696 0.085267 0.0060463 0.0046137];
-shear = [1.7094 0.55278 0.2054 0.085678];
-std_s = [0.2428 0.068223 0.012734 0.012291];
-subplot(121)
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.5$',...
-    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(5,:)); hold on;
-subplot(122)
-errorbar(eta,shear,std_s,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.5$',...
-    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(5,:)); hold on;
+% % 10,5 nuDGGK=0.01
+% eta    = [0.6,0.7,0.8 0.9];
+% gamma  = [2.4949 0.23368 0.057792 0.023572];
+% std_g  = [0.8696 0.085267 0.0060463 0.0046137];
+% shear = [1.7094 0.55278 0.2054 0.085678];
+% std_s = [0.2428 0.068223 0.012734 0.012291];
+% subplot(121)
+% errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.5$',...
+%     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(5,:)); hold on;
+% subplot(122)
+% errorbar(eta,shear,std_s,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.5$',...
+%     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(5,:)); hold on;
 
-% Mix len Ricci 2006
+kmax = [0.733 0.63 0.52 0.47];
+gmax = [0.32 0.22 0.11 0.027];
 subplot(121)
+% Mix len Ricci 2006
 semilogy([0.5  1.0],[10  1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
+% Bohm transport :
+btransp = gmax./kmax.^2;
+semilogy(.5:.1:.8,btransp','--','color',[0.4,0,0]+0.6, 'DisplayName','$\gamma/k^2$');
 set(gca, 'YScale', 'log'); grid on; legend('show')
 xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$') 
 subplot(122)
-grate = [0.2194 0.129 0.05084 0.01346];
-semilogy([0.6 0.7 0.8 0.9],grate,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
+semilogy([0.5:0.1:0.8],gmax,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
 set(gca, 'YScale', 'log'); grid on; legend('show')
 xlabel('$L_n/R$'); ylabel('$V_E$''') 
 
-%% nuDGGK = 0.1
+%% nuDGGK = 0.001
 figure; set(gcf, 'Position',  [100, 100, 1200, 350])
 % 6,3
-eta =   [0.6        0.7     0.8];
-gamma = [0.24321    0.085   0.0367];
-std_g  = [0.295      0.05    0.0023];
-gbmax  = [1.0        0.21    0.0367];
-gbmin  = [0.02       0.04    0.0367];
-% errorbar(eta,gamma,std_,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.1$',...
-%     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(6,:)); hold on;
+eta =   [0.6        0.7     0.8 0.9];
+gamma = [0.26 0.088 0.042 0.0156];
+std_g  = [0.4 0.06 0.003 0.0014];
+gbmax  = [1.2 0.33 0.042 0.0156];
+gbmin  = [0.015 0.035 0.003 0.0156];
+shear = [0.75 0.55 0.34 0.19];
+std_s  = [0.33 0.12 0.012];
+sbmax  = [1.4 0.796 0.34 0.19];
+sbmin  = [0.4 0.4 0.34 0.19];
+
 subplot(121)
-plot(eta,gamma,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
-    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(6,:)); hold on;
-plot(eta,gbmax,'-.','DisplayName',' ','Color',line_colors(6,:)); hold on;
-plot(eta,gbmin,'-.','DisplayName',' ','Color',line_colors(6,:)); hold on;
+plot(eta,gamma,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.001$',...
+    'LineWidth',2.0,'MarkerSize',5, 'Color', 'k'); hold on;
+subplot(122)
+plot(eta,shear,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.001$',...
+    'LineWidth',2.0,'MarkerSize',5, 'Color', 'k'); hold on;
 
 % 10,5
-eta =   [0.5        0.6        0.7      0.8       0.9];
-gamma = [12.2       0.19761 0.088 0.04253 0.026037];
-std_g   = [4.7      0.21328 0.065 0.009 0.00118];
-gbmax  = [12.2      0.8 0.25 0.04253 0.026037];
-gbmin  = [12.2      0.02 0.03 0.04253 0.026037];
-shear = [0.65361 0.46548 0.30645 0.21123 ];
-std_s  = [0.21288 0.10233 0.02886 0.0044664 ];
-sbmax  = [1 0.7 0.30645 0.21123 ];
-sbmin  = [0.4 0.35 0.30645 0.21123 ];
+% eta =   [0.5        0.6        0.7      0.8       0.9];
+% gamma = [12.2       0.19761 0.088 0.04253 0.026037];
+% std_g   = [4.7      0.21328 0.065 0.009 0.00118];
+% gbmax  = [12.2      0.8 0.25 0.04253 0.026037];
+% gbmin  = [12.2      0.02 0.03 0.04253 0.026037];
+% shear = [0.65361 0.46548 0.30645 0.21123 ];
+% std_s  = [0.21288 0.10233 0.02886 0.0044664 ];
+% sbmax  = [1 0.7 0.30645 0.21123 ];
+% sbmin  = [0.4 0.35 0.30645 0.21123 ];
 % errorbar(eta,gamma,std_,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
 %     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
-subplot(121)
-plot(eta,gamma,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
-    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
-plot(eta,gbmax,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
-plot(eta,gbmin,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
-subplot(122)
-plot(eta(2:end),shear,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
-    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
-plot(eta(2:end),sbmax,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
-plot(eta(2:end),sbmin,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
+% subplot(121)
+% plot(eta,gamma,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
+%     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
+% plot(eta,gbmax,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
+% plot(eta,gbmin,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
+% subplot(122)
+% plot(eta(2:end),shear,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
+%     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
+% plot(eta(2:end),sbmax,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
+% plot(eta(2:end),sbmin,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
 
-% Mix len Ricci 2006
+kmax = [0.27 0.68 0.52 0.31];
+gmax = [0.838 0.195 0.109 0.029];
 subplot(121)
+% Mix len Ricci 2006
 semilogy([0.5  1.0],[10  1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
+% Bohm transport :
+btransp = gmax./kmax.^2;
+semilogy(.6:.1:.9,btransp','--','color',[0.4,0,0]+0.6, 'DisplayName','$\gamma/k^2$');
 set(gca, 'YScale', 'log'); grid on; legend('show')
 xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$') 
 subplot(122)
-grate = [0.236 0.174 0.112 0.053];
-semilogy([0.6 0.7 0.8 0.9],grate,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
+semilogy([.6:.1:.9],gmax,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
 set(gca, 'YScale', 'log'); grid on; legend('show')
 xlabel('$L_n/R$'); ylabel('$V_E$''') 
-
-%% nuSGGK = 1.0
-figure
-% 6,3 nuDGGK=1.0
-eta = [0.5 0.6,0.7,0.8];
-gamma = [2.3 0.2215 0.0133 0.0032];
-std_g   = [3.1 0.22 0.0019 0.0006];
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=06,03$, $\nu_{SGGK}=1.0$',...
-    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(8,:)); hold on;
-
-
-% 10,5 nuDGGK = 1.0
-eta = [0.5 0.6,0.7,0.8];
-gamma = [10 0.319 0.009 0.0026];
-std_g   = [1.34 0.228 0.001 0.001];
-% errorbar(eta,gamma,err,'-.','DisplayName','$P,J=10,5$, $\nu_{SGGK}=1.0$',...
-%     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(1,:)); hold on;
-
-% Mix len Ricci 2006
-semilogy([0.5  1.0],[10  1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
-set(gca, 'YScale', 'log'); grid on; legend('show')
-xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$') 
+subplot(121)
+plot(eta,gbmax,'-.','DisplayName',' ','Color','k'); hold on;
+plot(eta,gbmin,'-.','DisplayName',' ','Color','k'); hold on;
+subplot(122)
+plot(eta,sbmax,'-.','DisplayName',' ','Color','k'); hold on;
+plot(eta,sbmin,'-.','DisplayName',' ','Color','k'); hold on;
-- 
GitLab