From cc9b95e7609c21068aa95aba11cc07a8c91aa1a1 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Wed, 8 Jun 2022 11:06:42 +0200
Subject: [PATCH] small script changes

---
 matlab/compute/mode_growth_meter.m            | 30 +++++++------
 .../plot_radial_transport_and_spacetime.m     | 42 ++++++++++---------
 wk/analysis_HeLaZ.m                           | 18 ++++----
 wk/analysis_gene.m                            | 16 +++----
 wk/header_3D_results.m                        |  3 +-
 5 files changed, 59 insertions(+), 50 deletions(-)

diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index b8c229d8..8c2bce8f 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -3,8 +3,12 @@ function [FIGURE] = mode_growth_meter(DATA,OPTIONS)
 NORMALIZED = OPTIONS.NORMALIZED;
 Nma   = OPTIONS.NMA; %Number moving average
 
-iz = OPTIONS.iz;
-[~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(1,:,iz,:))),2)));
+switch OPTIONS.iz
+    case 'avg'
+        field = squeeze(mean(DATA.PHI,3));
+    otherwise
+        field = squeeze(DATA.PHI(:,:,OPTIONS.iz,:));
+end
 
 FRAMES = zeros(size(OPTIONS.TIME));
 for i = 1:numel(OPTIONS.TIME)
@@ -13,6 +17,8 @@ end
 FRAMES = unique(FRAMES);
 t  = DATA.Ts3D(FRAMES);
 
+[~,ikzf] = max(mean(squeeze(abs(field(1,:,FRAMES))),2));
+
 FIGURE.fig = figure; set(gcf, 'Position',  [100 100 1200 700])
 FIGURE.FIGNAME = 'mode_growth_meter';
 for i = 1:2
@@ -20,27 +26,27 @@ for i = 1:2
 
     if MODES_SELECTOR == 2
         if NORMALIZED
-            plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES)))./max(abs(squeeze(x(1,ik,iz,FRAMES)))),Nma);
+            plt = @(x,ik) movmean(abs(squeeze(x(1,ik,FRAMES)))./max(abs(squeeze(x(1,ik,FRAMES)))),Nma);
         else
-            plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES))),Nma);
+            plt = @(x,ik) movmean(abs(squeeze(x(1,ik,FRAMES))),Nma);
         end
         kstr = 'k_x';
         k = DATA.kx;
         MODESTR = 'Zonal modes';
     elseif MODES_SELECTOR == 1
         if NORMALIZED
-            plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES)))./max(abs(squeeze(x(ik,1,iz,FRAMES)))),Nma);
+            plt = @(x,ik) movmean(abs(squeeze(x(ik,1,FRAMES)))./max(abs(squeeze(x(ik,1,FRAMES)))),Nma);
         else
-            plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES))),Nma);
+            plt = @(x,ik) movmean(abs(squeeze(x(ik,1,FRAMES))),Nma);
         end
         kstr = 'k_y';
         k = DATA.ky;
         MODESTR = 'NZ modes';
     elseif MODES_SELECTOR == 3
         if NORMALIZED
-            plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,iz,FRAMES)))./max(abs(squeeze(x(ik,ik,iz,FRAMES)))),Nma);
+            plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,FRAMES)))./max(abs(squeeze(x(ik,ik,FRAMES)))),Nma);
         else
-            plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,iz,FRAMES))),Nma);
+            plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,FRAMES))),Nma);
         end
         kstr = 'k_y=k_x';
         k = DATA.ky;
@@ -60,7 +66,7 @@ for i = 1:2
     amp   = MODES;
     i_=1;
     for ik = MODES
-        gr = polyfit(t,log(plt(DATA.PHI,ik)),1);
+        gr = polyfit(t,log(plt(field,ik)),1);
         gamma(i_) = gr(1);
         amp(i_)   = gr(2);
         i_=i_+1;
@@ -69,7 +75,7 @@ for i = 1:2
     %plot
     subplot(2,3,1+3*(i-1))
         [YY,XX] = meshgrid(t,fftshift(k));
-        pclr = pcolor(XX,YY,abs(plt(fftshift(DATA.PHI,MODES_SELECTOR),1:numel(k))));set(pclr, 'edgecolor','none');  hold on;
+        pclr = pcolor(XX,YY,abs(plt(fftshift(field,MODES_SELECTOR),1:numel(k))));set(pclr, 'edgecolor','none');  hold on;
         set(gca,'YDir','normal')
     %     xlim([t(1) t(end)]); %ylim([1e-5 1])
         xlabel(['$',kstr,'\rho_s$']); ylabel('$t c_s /\rho_s$');
@@ -78,11 +84,11 @@ for i = 1:2
     subplot(2,3,2+3*(i-1))
         mod2plot = [2:OPTIONS.NMODES+1];
         for i_ = mod2plot
-            semilogy(t,plt(DATA.PHI,MODES(i_))); hold on;
+            semilogy(t,plt(field,MODES(i_))); hold on;
     %         semilogy(t,exp(gamma(i_).*t+amp(i_)),'--k')
         end
         if MODES_SELECTOR == 2
-            semilogy(t,plt(DATA.PHI,ikzf),'--k');
+            semilogy(t,plt(field,ikzf),'--k');
         end
         xlim([t(1) t(end)]); %ylim([1e-5 1])
         xlabel('$t c_s /\rho_s$'); ylabel(['$|\phi_{',kstr,'}|$']);
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 93e0c5cb..d052fc15 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -3,12 +3,15 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
     tend = OPTIONS.TAVG_1; tstart = OPTIONS.TAVG_0;
     [~,its0D] = min(abs(DATA.Ts0D-tstart));
     [~,ite0D]   = min(abs(DATA.Ts0D-tend));
+    [~,its3D] = min(abs(DATA.Ts3D-tstart));
+    [~,ite3D]   = min(abs(DATA.Ts3D-tend));
     SCALE = DATA.scale;%(1/DATA.Nx/DATA.Ny)^2;
     Gx_infty_avg = mean(DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
     Gx_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
     Qx_infty_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE;
     Qx_infty_std = std (DATA.HFLUX_X(its0D:ite0D))*SCALE;
-    [~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(:,1,1,:))),2)));
+    f_avg_z      = squeeze(mean(DATA.PHI(:,:,:,:),3));
+    [~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3)));
     Ns3D = numel(DATA.Ts3D);
     [KX, KY] = meshgrid(DATA.kx, DATA.ky);
     
@@ -16,34 +19,33 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
 
     % Compute Gamma from ifft matlab
     Gx = zeros(DATA.Ny,DATA.Nx,numel(DATA.Ts3D));
-    for it = 1:numel(DATA.Ts3D)
-        for iz = 1:DATA.Nz
-            Gx(:,:,it)  = Gx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
-                          .*ifourier_GENE(DATA.DENS_I(:,:,iz,it));
-        end
-        Gx(:,:,it)  = Gx(:,:,it)/DATA.Nz;
-    end
+%     for it = 1:numel(DATA.Ts3D)
+%         for iz = 1:DATA.Nz
+%             Gx(:,:,it)  = Gx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
+%                           .*ifourier_GENE(DATA.DENS_I(:,:,iz,it));
+%         end
+%         Gx(:,:,it)  = Gx(:,:,it)/DATA.Nz;
+%     end
     Gx_t_mtlb = squeeze(mean(mean(Gx,1),2)); 
     % Compute Heat flux from ifft matlab
     Qx = zeros(DATA.Ny,DATA.Nx,numel(DATA.Ts3D));
-    for it = 1:numel(DATA.Ts3D)
-        for iz = 1:DATA.Nz
-            Qx(:,:,it)  = Qx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
-                          .*ifourier_GENE(DATA.TEMP_I(:,:,iz,it));
-        end
-        Qx(:,:,it)  = Qx(:,:,it)/DATA.Nz;
-    end
+%     for it = 1:numel(DATA.Ts3D)
+%         for iz = 1:DATA.Nz
+%             Qx(:,:,it)  = Qx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
+%                           .*ifourier_GENE(DATA.TEMP_I(:,:,iz,it));
+%         end
+%         Qx(:,:,it)  = Qx(:,:,it)/DATA.Nz;
+%     end
     Qx_t_mtlb = squeeze(mean(mean(Qx,1),2)); 
     % zonal vs nonzonal energies for phi(t)
 
     E_Zmode_SK       = zeros(1,Ns3D);
     E_NZmode_SK      = zeros(1,Ns3D);
     for it = 1:numel(DATA.Ts3D)
-        E_Zmode_SK(it)   = squeeze(DATA.ky(ikzf).^2.*abs(DATA.PHI(ikzf,1,1,it))^2);
-        E_NZmode_SK(it)  = squeeze(sum(sum(((1+KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2.*(KY~=0)))));
+        E_Zmode_SK(it)   = squeeze(DATA.ky(ikzf).^2.*abs(squeeze(f_avg_z(ikzf,1,it))).^2);
+        E_NZmode_SK(it)  = squeeze(sum(sum(((1+KX.^2+KY.^2).*abs(squeeze(f_avg_z(:,:,it))).^2.*(KY~=0)))));
     end
-    [~,its3D] = min(abs(DATA.Ts3D-tstart));
-    [~,ite3D]   = min(abs(DATA.Ts3D-tend));
+
 
 %% Figure    
 mvm = @(x) movmean(x,OPTIONS.NMVA);
@@ -83,7 +85,7 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
     toplot = process_field(DATA,OPTIONS);
     f2plot = toplot.FIELD;
     
-    clim = max(max(max(abs(plt(f2plot(:,:,:))))));
+    clim = max(max(max(abs(plt(f2plot(:,:,its3D:ite3D))))));
     subplot(313)
         [TY,TX] = meshgrid(DATA.x,DATA.Ts3D(toplot.FRAMES));
         pclr = pcolor(TX,TY,squeeze(plt(f2plot))'); 
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index 6fba72bd..4a06a7a1 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -10,7 +10,7 @@ system(['mkdir -p ',MISCDIR]);
 CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
 system(CMD);
 % Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 10; JOBNUMMAX = 20;
+JOBNUMMIN = 15; JOBNUMMAX = 16;
 data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
 data.localdir = LOCALDIR;
 data.FIGDIR   = LOCALDIR;
@@ -49,12 +49,12 @@ options.NAME      = '\phi';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'xz';
+options.PLAN      = 'xy';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
-options.COMP      = 'avg';
+options.COMP      = 9;
 % options.TIME      = dat.Ts5D;
-options.TIME      = 1250:1:1500;
+options.TIME      = [1800:1:2000];
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 create_film(data,options,'.gif')
@@ -63,7 +63,7 @@ end
 if 0
 %% 2D snapshots
 % Options
-options.INTERP    = 0;
+options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
 options.NAME      = '\phi';
@@ -72,11 +72,11 @@ options.NAME      = '\phi';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
-options.PLAN      = 'kxky';
+options.PLAN      = 'xy';
 % options.NAME      = 'f_i';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = [1200 1300 1400 1500];
+options.TIME      = [1600 1800 2000];
 data.a = data.EPS * 2e3;
 fig = photomaton(data,options);
 save_figure(data,fig)
@@ -165,10 +165,10 @@ if 0
 %% Mode evolution
 options.NORMALIZED = 0;
 options.K2PLOT = 1;
-options.TIME   = 1200:1300;
+options.TIME   = 1800:2000;
 options.NMA    = 1;
 options.NMODES = 15;
-options.iz     = 9;
+options.iz     = 'avg';
 fig = mode_growth_meter(data,options);
 save_figure(data,fig)
 end
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 281885de..86e064c5 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -16,7 +16,7 @@ options.TAVG_1   = gene_data.Ts3D(end); % Averaging times duration
 options.NMVA     = 1;              % Moving average for time traces
 options.ST_FIELD = '\phi';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x, Q_x)
 options.INTERP   = 1;
-gene_data.FIGDIR = folder;;
+gene_data.FIGDIR = folder;
 fig = plot_radial_transport_and_spacetime(gene_data,options);
 save_figure(gene_data,fig)
 end
@@ -32,11 +32,11 @@ options.NAME      = '\phi';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
-options.PLAN      = 'kxky';
+options.PLAN      = 'xy';
 % options.NAME      ='f_e';
 % options.PLAN      = 'sx';
-options.COMP      = 9;
-options.TIME      = [0 50 100 200 300];
+options.COMP      = 'avg';
+options.TIME      = [500];
 gene_data.a = data.EPS * 2000;
 fig = photomaton(gene_data,options);
 save_figure(gene_data,fig)
@@ -52,11 +52,11 @@ options.NAME      = '\phi';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'xz';
+options.PLAN      = 'xy';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = 000:300;
+options.TIME      = gene_data.Ts3D;
 gene_data.a = data.EPS * 2000;
 create_film(gene_data,options,'.gif')
 end
@@ -104,10 +104,10 @@ if 0
 %% Mode evolution
 options.NORMALIZED = 0;
 options.K2PLOT = 1;
-options.TIME   = 100:200;
+options.TIME   = 100:700;
 options.NMA    = 1;
 options.NMODES = 15;
-options.iz     = 9;
+options.iz     = 'avg';
 fig = mode_growth_meter(gene_data,options);
 save_figure(gene_data,fig)
 end
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 05a07248..9b790c18 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -5,7 +5,8 @@ helazdir = '/home/ahoffman/HeLaZ/';
 outfile ='';
 outfile ='';
 outfile ='';
-outfile ='shearless_cyclone/128x128x16xdmax_6_L_120_CBC_1.0';
+outfile ='shearless_cyclone/128x128x16x5x3_start';
+% outfile ='shearless_cyclone/128x128x16xdmax_6_L_120_CBC_1.0';
 % outfile ='shearless_cyclone/128x128x16xdmax_L_120_CBC_1.0';
 % outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8_to_1.0';
 % outfile ='quick_run/CLOS_1_64x64_5x3_L_120_kN_2.0_kT_0.5_nu_1e-01_SGGK';
-- 
GitLab