From 9dcecc98e468b0a5148a8e2f155db5723b401adf Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Tue, 5 Jul 2022 11:56:54 +0200
Subject: [PATCH] scripts updates

---
 matlab/plot/photomaton.m                      |  6 +--
 .../plot_radial_transport_and_spacetime.m     |  2 +-
 matlab/plot/spectrum_1D.m                     |  6 +--
 matlab/process_field.m                        | 20 ++++++++++
 wk/analysis_HeLaZ.m                           | 36 +++++++++---------
 wk/analysis_gene.m                            |  8 ++--
 wk/g2k3_transport_scaling.m                   | 38 +++++++++++++++++++
 wk/header_3D_results.m                        |  4 +-
 wk/quick_run.m                                | 24 +++++++-----
 9 files changed, 105 insertions(+), 39 deletions(-)
 create mode 100644 wk/g2k3_transport_scaling.m

diff --git a/matlab/plot/photomaton.m b/matlab/plot/photomaton.m
index fddaa5e6..c55123e0 100644
--- a/matlab/plot/photomaton.m
+++ b/matlab/plot/photomaton.m
@@ -27,9 +27,9 @@ FIGURE.fig = figure; set(gcf, 'Position',  toplot.DIMENSIONS.*[1 1 Ncols Nrows])
             if ~strcmp(OPTIONS.PLAN,'kxky')
                 caxis([-1,1]);
                 colormap(bluewhitered);
-            end
-            if OPTIONS.INTERP
-                shading interp; 
+                if OPTIONS.INTERP
+                    shading interp; 
+                end
             end
         else
             contour(toplot.X,toplot.Y,toplot.FIELD(:,:,i_)./scale,128);
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 5c21b786..a56fc536 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -68,7 +68,7 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
         plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_infty_avg, '-k',...
             'DisplayName',['$Q_x^{\infty} = $',num2str(Qx_infty_avg),'$\pm$',num2str(Qx_infty_std)]);
         ylabel('$Q_x$')  
-        ylim([0,2*abs(Qx_infty_avg)]); 
+%         ylim([0,2*abs(Qx_infty_avg)]); 
         xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
     grid on; set(gca,'xticklabel',[]); 
     title({DATA.param_title,...
diff --git a/matlab/plot/spectrum_1D.m b/matlab/plot/spectrum_1D.m
index f9ddf909..6c644301 100644
--- a/matlab/plot/spectrum_1D.m
+++ b/matlab/plot/spectrum_1D.m
@@ -68,8 +68,8 @@ if ~PLOT2D
         end
         plot(X,Y,'DisplayName','t-averaged')  
     otherwise
-    for it = 1:numel(toplot.FRAMES)
-        Gk    = compy(abs(toplot.FIELD(:,:,toplot.FRAMES(it))));
+    for it = 1:numel(frames)
+        Gk    = compx(abs(toplot.FIELD(:,:,it)));
         Gk    = squeeze(Gk);
         if options.NORM
             Gk    = Gk./max(max(abs(Gk)));
@@ -109,7 +109,7 @@ if ~PLOT2D
         plot(X,Y,'DisplayName','t-averaged')  
     otherwise
     for it = 1:numel(toplot.FRAMES)
-        Gk    = compx(abs(toplot.FIELD(:,:,it)));
+        Gk    = compy(abs(toplot.FIELD(:,:,it)));
         Gk    = squeeze(Gk);
         if options.NORM
             Gk    = Gk./max(max(abs(Gk)));
diff --git a/matlab/process_field.m b/matlab/process_field.m
index b6409f70..daa9f16e 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -176,6 +176,26 @@ switch OPTIONS.NAME
                 FIELD(:,:,it) = squeeze(compr(tmp));
             end                
         end
+    case 'N_e^{00}'
+        NAME = 'Ne00';
+        if COMPDIM == 3
+            for it = 1:numel(FRAMES)
+                tmp = squeeze(compr(DATA.Ne00(:,:,:,FRAMES(it))));
+                FIELD(:,:,it) = squeeze(process(tmp));
+            end
+        else
+            if REALP
+                tmp = zeros(Ny,Nx,Nz);
+            else
+                tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
+            end
+            for it = 1:numel(FRAMES)
+                for iz = 1:numel(DATA.z)
+                    tmp(:,:,iz) = squeeze(process(DATA.Ne00(:,:,iz,FRAMES(it))));
+                end
+                FIELD(:,:,it) = squeeze(compr(tmp));
+            end                
+        end
     case 'n_e'
         NAME = 'ne';
         if COMPDIM == 3
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index 5dc8efc3..b7d24578 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -11,7 +11,7 @@ system(['mkdir -p ',LOCALDIR]);
 CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
 system(CMD);
 % Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 20;
+JOBNUMMIN = 00; JOBNUMMAX = 02;
 data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
 data.localdir = LOCALDIR;
 data.FIGDIR   = LOCALDIR;
@@ -26,7 +26,7 @@ if 1
 options.TAVG_0   = 0.7*data.Ts3D(end); data.scale = (1/data.Nx/data.Ny)^2;
 options.TAVG_1   = data.Ts3D(end); % Averaging times duration
 options.NMVA     = 1;              % Moving average for time traces
-% options.ST_FIELD = '\Gamma_x';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+% options.ST_FIELD = '\Gamma_x';   % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
 options.ST_FIELD = '\phi';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
 options.INTERP   = 1;
 fig = plot_radial_transport_and_spacetime(data,options);
@@ -44,18 +44,18 @@ if 0
 % Options
 options.INTERP    = 1;
 options.POLARPLOT = 0;
-options.NAME      = '\phi';
+% options.NAME      = '\phi';
 % options.NAME      = 'N_i^{00}';
 % options.NAME      = 'v_y';
 % options.NAME      = 'n_i^{NZ}';
-% options.NAME      = '\Gamma_x';
+options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
 options.PLAN      = 'xy';
 % options.NAME      = 'f_i';
 % options.PLAN      = 'sx';
 options.COMP      = 1;
 % options.TIME      = data.Ts5D(end-30:end);
-options.TIME      =  data.Ts3D(1:2:end);
+options.TIME      =  data.Ts3D(450:end);
 % options.TIME      = [350:600];
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
@@ -65,20 +65,20 @@ end
 if 0
 %% 2D snapshots
 % Options
-options.INTERP    = 0;
+options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
-options.NAME      = '\phi';
+% options.NAME      = '\phi';
 % options.NAME      = 'n_i';
-% options.NAME      = 'N_i^{00}';
+% options.NAME      = 'N_e^{00}';
 % options.NAME      = 'T_i';
-% options.NAME      = '\Gamma_x';
+options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
 options.PLAN      = 'kxky';
 % options.NAME      'f_i';
 % options.PLAN      = 'sx';
-options.COMP      = 1;
-options.TIME      = [20 30 40 60 70];
+options.COMP      = 'avg';
+options.TIME      = [100 200 300 900 1300];
 data.a = data.EPS * 2e3;
 fig = photomaton(data,options);
 % save_figure(data,fig)
@@ -121,7 +121,7 @@ options.ST         = 1;
 options.PLOT_TYPE  = 'space-time';
 options.NORMALIZED = 1;
 options.JOBNUM     = 0;
-options.TIME       = [300:500];
+options.TIME       = [100:1200];
 options.specie     = 'i';
 options.compz      = 'avg';
 fig = show_moments_spectrum(data,options);
@@ -131,16 +131,16 @@ end
 
 if 0
 %% Time averaged spectrum
-options.TIME   = 650:800;
+options.TIME   = 1100:10:1400;
 options.NORM   =1;
-options.NAME   = '\phi';
-% options.NAME      = 'n_i';
-% options.NAME      ='\Gamma_x';
+% options.NAME   = '\phi';
+% options.NAME      = 'N_i^{00}';
+options.NAME      ='\Gamma_x';
 options.PLAN   = 'kxky';
 options.COMPZ  = 'avg';
 options.OK     = 0;
 options.COMPXY = 'avg';
-options.COMPT  = 'avg';
+options.COMPT  = 0;
 options.PLOT   = 'semilogy';
 fig = spectrum_1D(data,options);
 % save_figure(data,fig,'.png')
@@ -167,7 +167,7 @@ if 0
 %% Mode evolution
 options.NORMALIZED = 0;
 options.K2PLOT = 1;
-options.TIME   = 0:8000;
+options.TIME   = data.Ts3D;
 options.NMA    = 1;
 options.NMODES = 15;
 options.iz     = 'avg';
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 4f19bd12..4fc0cb77 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -6,10 +6,12 @@
 % folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.5/';
 % folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_1.0/';
 % folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_0.8/';
+% folder = '/misc/gene_results/HP_fig_2a_mu_1e-2/';
 % folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
 % folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
 % folder = '/misc/gene_results/LD_zpinch_1.6/';
-folder = '/misc/gene_results/ZP_16x8_kn_1.6_00/';
+% folder = '/misc/gene_results/ZP_16x8_kn_1.6/';
+folder = '/misc/gene_results/ZP_HP_kn_2.5/';
 gene_data = load_gene_data(folder);
 gene_data = invert_kxky_to_kykx_gene_results(gene_data);
 if 1
@@ -102,9 +104,9 @@ options.times   = 200:300;
 options.specie  = 'i';
 options.PLT_FCT = 'contour';
 options.folder  = folder;
-options.iz      = 9;
+options.iz      = 1;
 options.FIELD   = '<f_>';
-options.ONED    = 1;
+options.ONED    = 0;
 % options.FIELD   = 'Q_es';
 plot_fa_gene(options);
 end
diff --git a/wk/g2k3_transport_scaling.m b/wk/g2k3_transport_scaling.m
new file mode 100644
index 00000000..8e837768
--- /dev/null
+++ b/wk/g2k3_transport_scaling.m
@@ -0,0 +1,38 @@
+%% G ~ g^2/k^3 rule
+
+K_N =  [1.60 1.70 1.80 1.90 2.00 2.10 2.20 2.30 2.40 2.50];
+%DG
+g_DG   = [0.08 0.20 0.31 0.42 0.55 0.61 0.70 0.79 0.87 0.96];
+k_DG   = [0.21 0.31 0.37 0.42 0.47 0.47 0.47 0.52 0.52 0.47];
+G_DG   = g_DG.^2./k_DG.^3./K_N;
+%SG
+g_SG   = [0.08 0.19 0.30 0.40 0.49 0.59 0.68 0.76 0.85 0.93];
+k_SG   = [0.26 0.37 0.42 0.42 0.47 0.47 0.52 0.52 0.52 0.47];
+G_SG   = g_SG.^2./k_SG.^3./K_N;
+%LR
+g_LR   = [0.12 0.18 0.25 0.32 0.39 0.47 0.55 0.64 0.72 0.81];
+k_LR   = [0.31 0.31 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37];
+G_LR   = g_LR.^2./k_LR.^3./K_N;
+%LD
+g_LD   = [0.05 0.11 0.19 0.27 0.35 0.44 0.54 0.63 0.73 0.83];
+k_LD   = [0.26 0.26 0.31 0.31 0.31 0.31 0.31 0.31 0.31 0.31];
+G_LD   = g_LD.^2./k_LD.^3./K_N;
+% nu=0
+g_CL   = [0.22 0.27 0.32 0.37 0.44 0.51 0.59 0.67 0.76 0.85];
+k_CL   = [0.73 0.73 0.73 0.73 0.63 0.52 0.52 0.52 0.52 0.42];
+G_CL   = g_CL.^2./k_CL.^3./K_N;
+%% Plot
+
+figure;
+% plot(K_N,G_SG); hold on;
+% plot(K_N,G_DG); hold on;
+% plot(K_N,G_LR); hold on;
+% plot(K_N,G_LD); hold on;
+% plot(K_N,G_CL,'--k'); hold on;
+k = 0.3;
+plot(K_N,g_SG.^2/k.^3); hold on;
+plot(K_N,g_DG.^2/k.^3); hold on;
+plot(K_N,g_LR.^2/k.^3); hold on;
+plot(K_N,g_LD.^2/k.^3); hold on;
+plot(K_N,g_CL.^2/k.^3,'--k'); hold on;
+% plot(K_N,(G_DG+G_SG+G_LR+G_LD)/4);
\ No newline at end of file
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index f7de9f84..e703b678 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -15,6 +15,7 @@ helazdir = '/home/ahoffman/HeLaZ/';
 % outfile = 'shearless_cyclone/64x32x80x5x3_CBC_Npol_5_kine';
 % outfile = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine';
 % outfile = 'shearless_cyclone/64x32x160x5x3_CBC_Npol_10_kine';
+% outfile = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine';
 %% CBC
 % outfile ='shearless_cyclone/64x32x16x5x3_CBC_080';
 % outfile ='shearless_cyclone/64x32x16x5x3_CBC_scan/64x32x16x5x3_CBC_100';
@@ -26,5 +27,6 @@ helazdir = '/home/ahoffman/HeLaZ/';
 
 % outfile = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK';
 %% ZPINCH
-outfile ='Zpinch_rerun/Kn_2.5_200x48x5x3';
+% outfile ='Zpinch_rerun/Kn_2.5_200x48x5x3';
+outfile ='Zpinch_rerun/Kn_1.6_200x48x11x6';
 run analysis_HeLaZ
diff --git a/wk/quick_run.m b/wk/quick_run.m
index af5aee9b..55cc3c54 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -15,8 +15,8 @@ CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %% PHYSICAL PARAMETERS
 NU      = 0.0;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
-K_N     = 2.0;   % Density gradient drive
-K_T     = 0.5;   % Temperature '''
+K_N     = 2.5;   % Density gradient drive
+K_T     = 0.25*K_N;   % Temperature '''
 K_E     = 0.0;   % Electrostat '''
 % SIGMA_E = 0.05196152422706632;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
 SIGMA_E = 0.0233380;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
@@ -40,12 +40,12 @@ Q0      = 1.0;    % safety factor
 SHEAR   = 0.0;    % magnetic shear (Not implemented yet)
 EPS     = 0.0;    % inverse aspect ratio
 %% TIME PARMETERS
-TMAX    = 50;  % Maximal time unit
-DT      = 1e-2;   % Time step
-SPS0D   = 20;      % Sampling per time unit for 2D arrays
+TMAX    = 500;  % Maximal time unit
+DT      = 2e-2;   % Time step
+SPS0D   = 1;      % Sampling per time unit for 2D arrays
 SPS2D   = 0;      % Sampling per time unit for 2D arrays
-SPS3D   = 20;      % Sampling per time unit for 2D arrays
-SPS5D   = 1;    % Sampling per time unit for 5D arrays
+SPS3D   = 1;      % Sampling per time unit for 2D arrays
+SPS5D   = 1/2;    % Sampling per time unit for 5D arrays
 SPSCP   = 0;    % Sampling per time unit for checkpoints
 JOB2LOAD= -1;
 %% OPTIONS
@@ -53,8 +53,8 @@ SIMID   = 'dbg';  % Name of the simulation
 LINEARITY = 'linear';   % activate non-linearity (is cancelled if KXEQ0 = 1)
 % Collision operator
 % (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO      = 'DG';
-GKCO    = 0; % gyrokinetic operator
+CO      = 'SG';
+GKCO    = 1; % gyrokinetic operator
 ABCO    = 1; % interspecies collisions
 INIT_ZF = 0; ZF_AMP = 0.0;
 CLOS    = 0;   % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))s
@@ -110,6 +110,10 @@ if 1
 trange = [0.5 1]*data.Ts3D(end);
 nplots = 1;
 lg = compute_fluxtube_growth_rate(data,trange,nplots);
+[gmax,     kmax] = max(lg.g_ky(:,end));
+[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
+msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
+msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
 end
 
 if 0
@@ -172,4 +176,4 @@ figure
 plot(data.Ts3D(it0:it1), plt(data.PHI));
 xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
 title(sprintf('$k_x=$%2.2f, $k_y=0.00$',data.kx(ikx)))
-end
+end
\ No newline at end of file
-- 
GitLab