From 264438f6ab33d5b814f06dd7fe6e52be8503eac0 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Tue, 31 Jan 2023 17:22:08 +0100
Subject: [PATCH] Save the scripts

---
 matlab/compile_results.m                      | 23 ++++++-
 matlab/extract_fig_data.m                     |  2 +-
 matlab/load/load_3D_data.m                    | 32 +++-------
 matlab/load/load_gene_data.m                  |  2 +-
 matlab/plot/photomaton.m                      |  4 +-
 .../plot_radial_transport_and_spacetime.m     |  9 ++-
 matlab/plot/show_moments_spectrum.m           | 29 ++++++---
 matlab/plot/statistical_transport_averaging.m | 15 +++--
 wk/analysis_gene.m                            | 48 +++++++++-----
 wk/analysis_gyacomo.m                         | 38 ++++++-----
 wk/header_3D_results.m                        | 20 ++++--
 wk/p2_heatflux_salpha_convergence.m           | 64 +++++++++++++++----
 12 files changed, 190 insertions(+), 96 deletions(-)

diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index f6ccd44b..fccc7ee8 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -13,6 +13,7 @@ DATA.DT_EVOL  = []; %
 % FIELDS
 Nipj_    = []; Nepj_    = [];
 Ni00_    = []; Ne00_    = [];
+Nipjz_   = []; Nepjz_   = [];
 HFLUXI_   = [];
 HFLUXE_   = [];
 GGAMMAI_ = [];
@@ -154,6 +155,18 @@ while(CONTINUE)
             end
             [Ni00, Ts3D, ~] = load_3D_data(filename, 'Ni00');
              Ni00_ = cat(4,Ni00_,Ni00); clear Ni00
+            if KIN_E
+                try
+             Nepjz  = load_3D_data(filename, 'Nepjz');
+             Nepjz_ = cat(4,Nepjz_,Nepjz); clear Nepjz
+                catch
+                end
+            end
+%             try
+            [Nipjz, Ts3D, ~] = load_3D_data(filename, 'Nipjz');
+             Nipjz_ = cat(4,Nipjz_,Nipjz); clear Nipjz        
+%             catch
+%             end
         end
         if W_DENS
             if KIN_E
@@ -195,7 +208,9 @@ while(CONTINUE)
         Ts5D = [];
         if W_NAPJ
         [Nipj, Ts5D, ~] = load_5D_data(filename, 'moments_i');
+        tic
             Nipj_ = cat(6,Nipj_,Nipj); clear Nipj
+        toc
             if KIN_E
                 Nepj  = load_5D_data(filename, 'moments_e');
                 Nepj_ = cat(6,Nepj_,Nepj); clear Nepj
@@ -274,9 +289,11 @@ else
     DATA.scale = 1;%(1/Nx/Ny)^2;
     % Fields
     DATA.GGAMMA_RI = GGAMMAI_; DATA.PGAMMA_RI = PGAMMAI_; DATA.HFLUX_X = HFLUXI_;
-    DATA.Nipj = Nipj_; DATA.Ni00 = Ni00_; DATA.DENS_I = DENS_I_; DATA.TEMP_I = TEMP_I_;
+    DATA.Nipj = Nipj_; DATA.Ni00 = Ni00_; DATA.Nipjz = Nipjz_;
+    DATA.DENS_I = DENS_I_; DATA.TEMP_I = TEMP_I_;
     if(KIN_E)
-    DATA.Nepj = Nepj_; DATA.Ne00 = Ne00_; DATA.DENS_E = DENS_E_; DATA.TEMP_E = TEMP_E_;
+    DATA.Nepj = Nepj_; DATA.Ne00 = Ne00_; DATA.Nepjz = Nepjz_;
+    DATA.DENS_E = DENS_E_; DATA.TEMP_E = TEMP_E_;
     DATA.HFLUX_XE = HFLUXE_;
     end
     DATA.Ts5D = Ts5D_; DATA.Ts3D = Ts3D_; DATA.Ts0D = Ts0D_;
@@ -299,6 +316,8 @@ else
         num2str(DATA.Nx),'$, $(P,J)=(',num2str(DATA.PMAXI),',',...
         num2str(DATA.JMAXI),')$,',' $\mu_{hd}=$(',num2str(DATA.MUx),...
         ',',num2str(DATA.MUy),')'];
+    DATA.paramshort = [num2str(DATA.Pmaxi),'x',num2str(DATA.Jmaxi),'x',...
+        num2str(DATA.Nkx),'x',num2str(DATA.Nky),'x',num2str(DATA.Nz)];
     JOBNUM = LASTJOB;
 
     filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM);
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index 4be8de84..2e19518e 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -4,7 +4,7 @@
 % tw = [3000 4000];
 % tw = [4000 4500];
 % tw = [4500 5000];
-tw = [400 40000];
+tw = [450 650];
 
 fig = gcf;
 axObjs = fig.Children;
diff --git a/matlab/load/load_3D_data.m b/matlab/load/load_3D_data.m
index 37b54548..9c9df71d 100644
--- a/matlab/load/load_3D_data.m
+++ b/matlab/load/load_3D_data.m
@@ -1,31 +1,19 @@
 function [ data, time, dt ] = load_3D_data( filename, variablename )
-%LOAD_2D_DATA load a 2D variable stored in a hdf5 result file from HeLaZ
+%LOAD_3D_DATA load a 3D variable stored in a hdf5 result file from HeLaZ
     time     = h5read(filename,'/data/var3d/time');
-    kx       = h5read(filename,'/data/grid/coordkx');
-    ky       = h5read(filename,'/data/grid/coordky');
-    z        = h5read(filename,'/data/grid/coordz');
     dt    = h5readatt(filename,'/data/input','dt');
     cstart= h5readatt(filename,'/data/input','start_iframe3d'); 
     
-    data     = zeros(numel(ky),numel(kx),numel(z),numel(time));
-    [KX,KY]  = meshgrid(kx,ky);
+    % Find array size by loading the first output
+    tmp   = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+1,'%06d')]);
+    sz    = size(tmp.real);
+    % add time dimension
+    sz    = [sz numel(time)];
+    data     = zeros(sz);
     
-    switch variablename
-    case 'vEx'
-         for it = 1:numel(time)
-            tmp         = h5read(filename,['/data/var3d/phi/', num2str(cstart+it,'%06d')]);
-            data(:,:,:,it) = 1i.*KY.*(tmp.real + 1i * tmp.imaginary);
-         end
-    case 'vEy'
-         for it = 1:numel(time)
-            tmp         = h5read(filename,['/data/var3d/phi/', num2str(cstart+it,'%06d')]);
-            data(:,:,:,it) = -1i.*KX.*(tmp.real + 1i * tmp.imaginary);
-        end           
-    otherwise
-        for it = 1:numel(time)
-            tmp         = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
-            data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
-        end
+    for it = 1:numel(time)
+        tmp         = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
+        data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
     end
 
 end
diff --git a/matlab/load/load_gene_data.m b/matlab/load/load_gene_data.m
index 63554d8a..2fdff459 100644
--- a/matlab/load/load_gene_data.m
+++ b/matlab/load/load_gene_data.m
@@ -1,6 +1,6 @@
 function [ DATA ] = load_gene_data( folder )
 %to load gene data as for HeLaZ results
-namelist      = read_namelist([folder,'parameters.dat']);
+namelist      = read_namelist([folder,'parameters']);
 DATA.namelist = namelist;
 DATA.folder   = folder;
 %% Grid
diff --git a/matlab/plot/photomaton.m b/matlab/plot/photomaton.m
index e60fd8cc..2e7ccdcf 100644
--- a/matlab/plot/photomaton.m
+++ b/matlab/plot/photomaton.m
@@ -15,8 +15,8 @@ Ncols  = ceil(Nframes/Nrows);
 %
 TNAME = [];
 FIGURE.fig = figure; set(gcf, 'Position',  toplot.DIMENSIONS.*[1 1 Ncols Nrows])
-    frame_max = max(max(max(abs(toplot.FIELD(:,:,:)))));
     for i_ = 1:numel(FRAMES)
+    frame_max = max(max(max(abs(toplot.FIELD(:,:,i_)))));
     subplot(Nrows,Ncols,i_); TNAME = [TNAME,'_',sprintf('%.0f',DATA.Ts3D(FRAMES(i_)))];
     if OPTIONS.NORMALIZE
         scale = frame_max; % Scaling to normalize
@@ -43,7 +43,7 @@ FIGURE.fig = figure; set(gcf, 'Position',  toplot.DIMENSIONS.*[1 1 Ncols Nrows])
         end
         xlabel(toplot.XNAME); ylabel(toplot.YNAME);
 %         if i_ > 1; set(gca,'ytick',[]); end; 
-        title([sprintf('$t c_s/R=%.0f$',tshot),', max = ',sprintf('%.1e',scale)]);
+        title([sprintf('$t c_s/R=%.0f$',tshot),', max = ',sprintf('%.1e',frame_max)]);
 
     end
     legend(['$',toplot.FIELDNAME,'$']);
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index e37e1073..884a919e 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -61,11 +61,16 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS,CODE)
     end
 
 %% Figure    
+clr_ = lines(20);
 mvm = @(x) movmean(x,OPTIONS.NMVA);
     FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; %set(gcf, 'Position',  [500, 1000, 1000, 600])
     FIGURE.ax1 = subplot(3,1,1,'parent',FIGURE.fig);
-        plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
-        plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
+        plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'--',...
+            'color',clr_((DATA.Pmaxi-1)/2-1,:),...
+            'DisplayName',['$\Gamma_x$ ',DATA.paramshort]); hold on;
+        plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X*SCALE),'-',...
+            'color',clr_((DATA.Pmaxi-1)/2-1,:),...
+            'DisplayName',['$Q_x$ ',DATA.paramshort]); hold on;
         ylabel('Transport')  
         if(~isnan(Qx_infty_avg))
         plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_infty_avg, '-k',...
diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m
index 2fc32fa4..3f000729 100644
--- a/matlab/plot/show_moments_spectrum.m
+++ b/matlab/plot/show_moments_spectrum.m
@@ -1,16 +1,25 @@
 function [ FIGURE ] = show_moments_spectrum( DATA, OPTIONS )
 
-
 Pi = DATA.Pi;
 Ji = DATA.Ji;
-Nipj = sum(sum(sum(abs(DATA.Nipj),3),4),5);
+if (isempty(DATA.Nipjz))
+    Time_ = DATA.Ts3D;
+    Nipj = sum(abs(DATA.Nipjz),3);
+else
+    Time_ = DATA.Ts5D;
+    Nipj = sum(sum(sum(abs(DATA.Nipj),3),4),5);
+end
 Nipj = squeeze(Nipj);
 
 if DATA.KIN_E
 Pe = DATA.Pe;
 Je = DATA.Je;
-Nepj = sum(sum(sum(abs(DATA.Nepj),3),4),5);
-Nepj = squeeze(Nepj);
+    if ~(isempty(DATA.Nepjz))
+        Nepj = sum(abs(DATA.Nepjz),3);
+    else
+        Nepj = sum(sum(sum(abs(DATA.Nepj),3),4),5);
+    end
+    Nepj = squeeze(Nepj);
 end
 
 FIGURE.fig = figure; FIGURE.FIGNAME = ['mom_spectrum_',DATA.PARAMS];
@@ -19,8 +28,8 @@ set(gcf, 'Position',  [100 50 1000 400])
 if ~OPTIONS.ST
     t0 = OPTIONS.TIME(1);
     t1 = OPTIONS.TIME(end);
-    [~,it0] = min(abs(t0-DATA.Ts5D));
-    [~,it1] = min(abs(t1-DATA.Ts5D));
+    [~,it0] = min(abs(t0-Time_));
+    [~,it1] = min(abs(t1-Time_));
     Nipj = mean(Nipj(:,:,it0:it1),3);
     if DATA.KIN_E
     Nepj = mean(Nepj(:,:,it0:it1),3);
@@ -87,7 +96,7 @@ else
     % weights to average
     weights = zeros(size(p2ji));
     % space time of moments amplitude wrt to p+2j degree
-    Ni_ST = zeros(numel(p2ji),numel(DATA.Ts5D));
+    Ni_ST = zeros(numel(p2ji),numel(Time_));
     % fill the st diagramm by averaging every p+2j=n moments together
     for ip = 1:numel(Pi)
         for ij = 1:numel(Ji)
@@ -109,7 +118,7 @@ else
     % weights to average
     weights = zeros(size(p2je));
     % space time of moments amplitude wrt to p+2j degree
-    Ne_ST = zeros(numel(p2je),numel(DATA.Ts5D));
+    Ne_ST = zeros(numel(p2je),numel(Time_));
     % fill the st diagramm by averaging every p+2j=n moments together
     for ip = 1:numel(Pe)
         for ij = 1:numel(Je)
@@ -134,7 +143,7 @@ else
     if DATA.KIN_E
     subplot(2,1,1)
     end
-        imagesc(DATA.Ts5D,p2ji,plt(Ni_ST,1:numel(p2ji))); 
+        imagesc(Time_,p2ji,plt(Ni_ST,1:numel(p2ji))); 
         set(gca,'YDir','normal')        
 %         pclr = pcolor(XX,YY,plt(Ni_ST));
 %         set(pclr, 'edgecolor','none'); hold on;
@@ -142,7 +151,7 @@ else
     title('$\langle\sum_k |N_i^{pj}|\rangle_{p+2j=const}$')
     if DATA.KIN_E
     subplot(2,1,2)
-        imagesc(DATA.Ts5D,p2je,plt(Ne_ST,1:numel(p2ji))); 
+        imagesc(Time_,p2je,plt(Ne_ST,1:numel(p2ji))); 
         set(gca,'YDir','normal')
 %         pclr = pcolor(XX,YY,plt(Ne_ST)); 
 %         set(pclr, 'edgecolor','none'); hold on;
diff --git a/matlab/plot/statistical_transport_averaging.m b/matlab/plot/statistical_transport_averaging.m
index ffc0395e..af463893 100644
--- a/matlab/plot/statistical_transport_averaging.m
+++ b/matlab/plot/statistical_transport_averaging.m
@@ -1,5 +1,5 @@
 
-function [ fig, Gx_infty_avg, Gx_infty_std ] = statistical_transport_averaging( data, options )
+function [ fig, res ] = statistical_transport_averaging( data, options )
 scale = data.scale;
 Trange  = options.T;
 [~,it0] = min(abs(Trange(1)-data.Ts0D)); 
@@ -21,7 +21,7 @@ dt_const = numel(unique(round(diff(data.Ts0D(it0:it1))*100)))==1;
         transp_seg_std(Np) = std(gamma(1:Np));
     end
 
-    time_seg = (data.Ts0D(it0:it1)-data.Ts0D(it0)); 
+    time_seg = (data.Ts0D(it0:it1)); 
     
     fig = 0;
 if options.NPLOTS > 0
@@ -43,8 +43,13 @@ if options.NPLOTS > 0
     xlabel('Averaging time'); ylabel('$\langle Q_x\rangle_{\tau}$');
     legend(['$Q_x^\infty=$',sprintf('%2.2e',transp_seg_avg(end))])
 end   
-Gx_infty_avg = mean(gamma);
-Gx_infty_std = std (gamma);
-disp(['G_x=',sprintf('%2.2e',Gx_infty_avg),'+-',sprintf('%2.2e',Gx_infty_std)]);
+res.time_seg = time_seg;
+res.Qx_t     = 
+res.Gx_avg = mean(gamma);
+res.Gx_std = std (gamma);
+disp(['G_x=',sprintf('%2.2e',res.Gx_avg),'+-',sprintf('%2.2e',res.Gx_std)]);
+res.Qx_avg = mean(Qx);
+res.Qx_std = std (Qx);
+disp(['G_x=',sprintf('%2.2e',res.Qx_avg),'+-',sprintf('%2.2e',res.Qx_std)]);
 end
 
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 0c3c014d..11c0b16a 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -29,7 +29,7 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
 % folder = '/misc/gene_results/CBC/128x64x16x24x12/';
 % folder = '/misc/gene_results/CBC/196x96x20x32x16_01/';
 % folder = '/misc/gene_results/CBC/256x96x24x32x12/';
-folder = '/misc/gene_results/CBC/128x64x16x6x4/';
+% folder = '/misc/gene_results/CBC/128x64x16x6x4/';
 % folder = '/misc/gene_results/CBC/KT_5.3_128x64x16x24x12_01/';
 % folder = '/misc/gene_results/CBC/KT_4.5_128x64x16x24x12_01/';
 % folder = '/misc/gene_results/CBC/KT_9_128x64x16x24x12/';
@@ -46,10 +46,17 @@ folder = '/misc/gene_results/CBC/128x64x16x6x4/';
 
 %Paper 2
 % folder = '/misc/gene_results/CBC/KT_6.96_64x32x32x24x12_Nexc_5/';
+% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x8x4_Nexc_5_00/';
+% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x16x8_Nexc_5_00/';
 % folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x32x16_Nexc_5_00/';
 % folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x32x16_Nexc_5_01/';
+
 % folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x32x16_Nexc_5_00/';
 % folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x32x16_Nexc_5_01/';
+% folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x16x8_Nexc_5/';
+% folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x8x4_Nexc_5/';
+folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x8x4_Nexc_5_smallvbox/';
+% folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x16x8_Nexc_5_largexbox/';
 
 gene_data = load_gene_data(folder);
 gene_data.FIGDIR = folder;
@@ -82,24 +89,31 @@ fig = statistical_transport_averaging(gene_data,options);
 end
 
 if 0
-%% 2D snapshots
+%% fields snapshots
 % Options
-options.INTERP    = 1;
+options.INTERP    = 0;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 0;
+options.NORMALIZE = 0;
+% options.NAME      = '\phi';
+% options.NAME      = '\psi';
+% options.NAME      = '\omega_z';
+options.NAME      = 'n_i';
+% options.NAME      = 'n_i-n_e';
+% options.NAME      = '\phi^{NZ}';
+% options.NAME      = 'N_i^{00}';
+% options.NAME      = 'N_i^{00}-N_e^{00}';
+% options.NAME      = 's_{Ex}';
 % options.NAME      = 'Q_x';
-options.NAME      = '\phi';
-% options.NAME      = 'n_i';
-% options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
-options.PLAN      = 'kxz';
-% options.NAME      ='f_e';
-% options.PLAN      = 'sx';
-options.COMP      = 'avg';
-options.TIME      = [40 55 70];
-gene_data.a = data.EPS * 2000;
+options.PLAN      = 'xy';
+options.COMP      = 1;
+options.TIME      = [0];
+options.RESOLUTION = 256;
+
+data.a = data.EPS * 2e3;
 fig = photomaton(gene_data,options);
-save_figure(gene_data,fig,'.png')
+% save_figure(data,fig)
 end
 
 if 0
@@ -111,7 +125,7 @@ options.NAME      = '\phi';
 % options.NAME      = 'v_{Ey}';
 % options.NAME      = 'G_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'xz';
+options.PLAN      = 'xy';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
@@ -168,11 +182,11 @@ end
 
 if 0
 %% Time averaged spectrum
-options.TIME   = [100 500];
+options.TIME   = [1];
 options.NORM   =1;
-% options.NAME   = '\phi';
+options.NAME   = '\phi';
 % options.NAME      = 'n_i';
-options.NAME   ='\Gamma_x';
+% options.NAME   ='\Gamma_x';
 options.PLAN   = 'kxky';
 options.COMPZ  = 'avg';
 options.OK     = 0;
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index a7f93e8a..000f5022 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -11,12 +11,13 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
 %% Load the results
 LOCALDIR  = [gyacomodir,resdir,'/'];
 MISCDIR   = ['/misc/gyacomo_outputs/',resdir,'/']; %For long term storage
+DATADIR   = [PARTITION,resdir,'/'];
 system(['mkdir -p ',MISCDIR]);
 system(['mkdir -p ',LOCALDIR]);
 % CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
 % system(CMD);
 % Load outputs from jobnummin up to jobnummax
-data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
+data = compile_results(DATADIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
 data.localdir = LOCALDIR;
 data.FIGDIR   = LOCALDIR;
 data.folder   = LOCALDIR;
@@ -46,17 +47,22 @@ end
 
 if 0
 %% statistical transport averaging
-gavg =[]; gstd = [];
+Gavg =[]; Gstd = [];
+Qavg =[]; Qstd = [];
+figure; hold on;
+plot(data.Ts0D,data.Qx);
 for i_ = 1:2:numel(data.TJOB_SE) 
 % i_ = 3; 
 disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))])
 disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))])
-options.T = [data.TJOB_SE(i_) data.TJOB_SE(i_+1)];
+options.T = [data.TJOB_SE(i_)*1.2 data.TJOB_SE(i_+1)];
 options.NPLOTS = 0;
-[fig, gavg_, gstd_] = statistical_transport_averaging(data,options);
-gavg = [gavg gavg_]; gstd = [gstd gstd_];
+[fig, res] = statistical_transport_averaging(data,options);
+Gavg = [Gavg res.Gx_avg]; Gstd = [Gstd res.Gx_std];
+Qavg = [Qavg res.Qx_avg]; Qstd = [Qstd res.Qx_std];
 end
-disp(gavg); disp(gstd);
+% disp(Gavg); disp(Gstd);
+disp(Qavg); disp(Qstd);
 end
 if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -66,18 +72,18 @@ options.POLARPLOT = 0;
 options.NAME      = '\phi';
 % options.NAME      = '\omega_z';
 % options.NAME     = 'N_i^{00}';
-% options.NAME      = 'v_x';
+% options.NAME      = 's_{Ey}';
 % options.NAME      = 'n_i^{NZ}';
-% options.NAME      = '\Gamma_x';
+% options.NAME      = 'G_x';
 % options.NAME      = 'n_i';
 % options.NAME      = 'n_i-n_e';
 options.PLAN      = 'xy';
-options.NAME      = 'f_i';
-options.PLAN      = 'sx';
+% options.NAME      = 'f_i';
+% options.PLAN      = 'sx';
 options.COMP      = 'avg';
 % options.TIME      = data.Ts5D(end-30:end);
-% options.TIME      =  data.Ts3D;
-options.TIME      = [0:10000];
+options.TIME      =  data.Ts3D;
+% options.TIME      = [0:10000];
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 options.RESOLUTION = 256;
@@ -94,7 +100,7 @@ options.NORMALIZE = 0;
 options.NAME      = '\phi';
 % options.NAME      = '\psi';
 % options.NAME      = '\omega_z';
-% options.NAME      = 'n_e';
+% options.NAME     = 'n_i';
 % options.NAME      = 'n_i-n_e';
 % options.NAME      = '\phi^{NZ}';
 % options.NAME      = 'N_i^{00}';
@@ -103,8 +109,8 @@ options.NAME      = '\phi';
 % options.NAME      = 'Q_x';
 % options.NAME      = 'k^2n_e';
 options.PLAN      = 'xy';
-options.COMP      = 'avg';
-options.TIME      = [500];
+options.COMP      = 1;
+options.TIME      = [29.5 30 30.5];
 options.RESOLUTION = 256;
 
 data.a = data.EPS * 2e3;
@@ -150,7 +156,7 @@ if 0
 %% Hermite-Laguerre spectrum
 % options.TIME = 'avg';
 options.P2J        = 0;
-options.ST         = 1;
+options.ST         = 0;
 options.PLOT_TYPE  = 'space-time';
 options.NORMALIZED = 0;
 options.JOBNUM     = 0;
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 9a860e05..9d8b87fe 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -1,9 +1,8 @@
 % Directory of the code "mypathtoHeLaZ/HeLaZ/"
 gyacomodir = '/home/ahoffman/gyacomo/';
-% Directory of the simulation (from results)
-% if 1% Local results
-% resdir ='volcokas/64x32x16x5x3_kin_e_npol_1';
-
+% Partition of the computer where the data have to be searched
+PARTITION  = '/misc/gyacomo_outputs/';
+% PARTITION  = gyacomodir;
 %% Dimits
 % resdir ='shearless_cyclone/128x64x16x5x3_Dim_90';
 % resdir ='shearless_cyclone/128x64x16x9x5_Dim_scan/128x64x16x9x5_Dim_60';
@@ -79,6 +78,7 @@ gyacomodir = '/home/ahoffman/gyacomo/';
 % resdir = 'GCM_CBC/daint/Miller_GX_comparison';
 % resdir = 'GCM_CBC/daint/Salpha_GX_comparison';
 %% Paper 2 simulations
+% convergence CBC and Dimits regime
 % resdir = 'paper_2_nonlinear/kT_6.96/CBC_7x4x128x64x24';
 % resdir = 'paper_2_nonlinear/kT_6.96/CBC_3x2x128x64x24_Nexc_5';
 % resdir = 'paper_2_nonlinear/kT_6.96/CBC_5x3x128x64x24_Nexc_5';
@@ -90,9 +90,17 @@ gyacomodir = '/home/ahoffman/gyacomo/';
 % resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_7x4x128x64x24_Nexc_5';
 % resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_9x5x128x64x24_Nexc_5';
 % resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_11x6x128x64x24_Nexc_5';
+% resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_13x7x128x64x24_Nexc_5';
 % resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_15x8x128x64x24_Nexc_5';
-resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_17x9x128x64x24_Nexc_5';
+% resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_17x9x128x64x24_Nexc_5';
+% Scan in kT
+% resdir = 'paper_2_nonlinear/kT_scan_DGGK_0.05/9x5x128x64x24';
+
+% resdir = 'paper_2_nonlinear/CBC_rerun/rerun_CBC_3x2x128x64x16';
+% resdir = 'dev/init_ppj';
+% resdir = 'dev/hatB_NL';
+resdir = 'dev/CBC_wave_study';
 
 resdir = ['results/',resdir];
-JOBNUMMIN = 00; JOBNUMMAX = 10;
+JOBNUMMIN = 04; JOBNUMMAX = 10;
 run analysis_gyacomo
diff --git a/wk/p2_heatflux_salpha_convergence.m b/wk/p2_heatflux_salpha_convergence.m
index 2d850670..49773c32 100644
--- a/wk/p2_heatflux_salpha_convergence.m
+++ b/wk/p2_heatflux_salpha_convergence.m
@@ -1,19 +1,59 @@
+%% Heat flux convergence for kt=6.96 and 5.3
 figure
-%% KT 6.96, nuDGDK = 0.05, 128x64x24, Nexc 5
+title('s-$\alpha$ turb. heat flux conv.');
+% KT 6.96, nuDGDK = 0.05, 128x64x16, Nexc 1
+P   = [2     4     12];
+Qx  = [50.00 46.00 41.00];
+std_= [6.600 2.300 6.600];
+    errorbar(P,Qx,std_/2,'o-.r',...
+    'LineWidth',2.0,'MarkerSize',8,...
+    'DisplayName','KT 6.9, nuDGDK 0.05'); hold on
+xlabel('$P$, $J=P/2$'); ylabel('$Q_x$');
+% KT 6.96, nuDGDK = 0.05, 128x64x24, Nexc 5
 P   = [4     6     8     10   ];
 Qx  = [67.62 67.50 59.21 64.17];
-std = [15.42 20.32 17.25 16.05];
-    errorbar(P,Qx,std/2,'s-',...
-    'LineWidth',2.0,...
-    'DisplayName','KT 6.96, nuDGDK 0.05'); hold on
+std_= [15.42 20.32 17.25 16.05];
+    errorbar(P,Qx,std_/2,'o-r',...
+    'LineWidth',2.0,'MarkerSize',8,...
+    'DisplayName','KT 6.9, nuDGDK 0.05'); hold on
 xlabel('$P$, $J=P/2$'); ylabel('$Q_x$');
-%% KT 5.3, nuDGDK = 0.05, 128x64x24, Nexc 5
-P   = [4     6     8     10   ];
-Qx  = [44.10 21.61 16.04 0.489];
-std = [10.61 6.952 4.166 0.061];
-    errorbar(P,Qx,std/2,'s-',...
-    'LineWidth',2.0,...
+% KT 5.3, nuDGDK = 0.05, 128x64x24, Nexc 5
+P   = [4     6     8     10    12   ];
+Qx  = [44.10 21.61 16.04 0.558 0.901];
+std_= [10.61 6.952 4.166 0.025 0.122];
+errorbar(P,Qx,std_/2,'o-b',...
+    'LineWidth',2.0,'MarkerSize',8,...
     'DisplayName','KT 5.3, nuDGDK 0.05');
 xlabel('$P$, $J=P/2$'); ylabel('$Q_x$');
 
-title('GYAC, turb. heat flux conv.');
\ No newline at end of file
+% GENE RESULTS 6.96 128x64x24, Nexc 5
+Nvp = [32    16    8     8];
+Qx  = [34.53 37.96 1.948 13.0];
+std_= [7.830 6.048 0.629 3.09];
+errorbar(Nvp,Qx,std_/2,'s--r',...
+    'LineWidth',2.0,'MarkerSize',8,...
+    'DisplayName','KT 6.9 GENE');
+xlabel('$P=N_{v\parallel}$, $J=P/2=N_\mu$'); ylabel('$Q_x$');
+
+% GENE RESULTS 5.3 128x64x24, Nexc 5
+% Comment: the result with nvp = 8 is not trusworthy as GENE does not have
+% a linear instability with this resolution.. It is then hard to draw
+% any conclusion from it.
+Nvp = [32    16    8];
+Qx  = [0.284 0.000 0.370];
+std_= [0.177 0.000 0.140];
+errorbar(Nvp,Qx,std_/2,'s--b',...
+    'LineWidth',2.0,'MarkerSize',8,...
+    'DisplayName','KT 5.3 GENE');
+xlabel('$P=N_{v\parallel}$, $J=P/2=N_\mu$'); ylabel('$Q_x$');
+
+%% KT scans 9x5x128x64x24
+clrs = lines(10);
+figure
+kT_ = [6.96      6.3       5.8        5.3         4.8];
+Qx_ = [59.21     40.1740   31.3627    16.04       1.0900];
+std_= [6.7880    7.5513    9.7815     4.166       0.4995];
+errorbar(kT_,Qx_,std_/2,'s--','color',clrs(4,:),...
+    'LineWidth',2.0,'MarkerSize',8,...
+    'DisplayName','(8,4)');
+xlabel('$K_T$'); ylabel('$Q_x$');
-- 
GitLab