From a9696f933d112bd1dd2945c1521b5a97a0f5b6fe Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Mon, 18 Jul 2022 15:34:29 +0200
Subject: [PATCH] save latest scripts versions

---
 matlab/compute/compute_fluxtube_growth_rate.m |  6 +-
 matlab/evolve_tracers.m                       | 37 ++++++++-----
 matlab/load/load_gene_data.m                  | 18 +++---
 matlab/load/quick_gene_load.m                 |  2 +-
 matlab/plot/plot_ballooning.m                 | 55 ++++++++++++-------
 matlab/plot/spectrum_1D.m                     |  3 +
 wk/analysis_HeLaZ.m                           | 16 +++---
 wk/analysis_gene.m                            | 15 ++---
 wk/g2k3_transport_scaling.m                   | 10 ++--
 wk/header_3D_results.m                        |  7 ++-
 wk/quick_run.m                                | 49 +++++++++--------
 11 files changed, 125 insertions(+), 93 deletions(-)

diff --git a/matlab/compute/compute_fluxtube_growth_rate.m b/matlab/compute/compute_fluxtube_growth_rate.m
index 9d93e421..47ba061c 100644
--- a/matlab/compute/compute_fluxtube_growth_rate.m
+++ b/matlab/compute/compute_fluxtube_growth_rate.m
@@ -42,9 +42,9 @@ linear_gr.ky     = kys;
 if PLOT >0
        figure
        plot(linear_gr.ky,linear_gr.g_ky(:,end),'-o','DisplayName','$Re(\omega_{k_y})$'); hold on;
-       plot(linear_gr.ky,linear_gr.w_ky(:,end),'-o','DisplayName','$Im(\omega_{k_y})$'); hold on;
-       plot(linear_gr.ky,linear_gr.ce  (:,end),'-o','DisplayName','$\epsilon$'); hold on;
-       xlim([min(linear_gr.ky) max(linear_gr.ky)]);
+       plot(linear_gr.ky,linear_gr.w_ky(:,end),'--*','DisplayName','$Im(\omega_{k_y})$'); hold on;
+       plot(linear_gr.ky,linear_gr.ce  (:,end),'x','DisplayName','$\epsilon$'); hold on;
+%        xlim([min(linear_gr.ky) max(linear_gr.ky)]);
        xlabel('$k_y$');
        legend('show');
        title(DATA.param_title);
diff --git a/matlab/evolve_tracers.m b/matlab/evolve_tracers.m
index c2616287..952a87ad 100644
--- a/matlab/evolve_tracers.m
+++ b/matlab/evolve_tracers.m
@@ -1,18 +1,19 @@
 % Options
 SHOW_FILM = 0;
+field2plot  ='Ni00';
 INIT     = 'lin';   % lin (for a line)/ round (for a small round)/ gauss for random
-U_TIME   = 700;     % >0 for frozen velocity at a given time, -1 for evolving field
+U_TIME   = 200;     % >0 for frozen velocity at a given time, -1 for evolving field
 Evolve_U = 1;       % 0 for frozen velocity at a given time, 1 for evolving field
-Tfin     = 10;
-dt_      = 0.0005;
+Tfin     = 200;
+dt_      = 0.1;
 Nstep    = ceil(Tfin/dt_);
 % Init tracers
-Np      = 200; %number of tracers
+Np      = 20; %number of tracers
 % color = tcolors;
 color = jet(Np);
 tcolors = distinguishable_colors(Np); %Their colors
 
-Na = 1000000; %length of trace
+Na = 1/dt_; %length of trace
 
 Traj_x = zeros(Np,Nstep);
 Traj_y = zeros(Np,Nstep);
@@ -60,9 +61,9 @@ for iz = 1:data.Nz
 %     Ux(:,:,iz) = real(ifft2( 1i*KY.*(data.PHI(:,:,iz,itu_)),data.Nx,data.Ny));
 %     Uy(:,:,iz) = real(ifft2(-1i*KX.*(data.PHI(:,:,iz,itu_)),data.Nx,data.Ny));
 %     ni(:,:,iz) = real(ifft2(data.DENS_I(:,:,iz,itu_),data.Nx,data.Ny));
-    Ux(:,:,iz) = ifourier_GENE( 1i*KY.*(data.PHI(:,:,iz,itu_)))'*sqrt(data.scale);
-    Uy(:,:,iz) = ifourier_GENE(-1i*KX.*(data.PHI(:,:,iz,itu_)))'*sqrt(data.scale);
-    ni(:,:,iz) = ifourier_GENE(data.DENS_I(:,:,iz,itu_))'*sqrt(data.scale);
+    Ux(:,:,iz) = ifourier_GENE( 1i*KY.*(data.PHI(:,:,iz,itu_)))';
+    Uy(:,:,iz) = ifourier_GENE(-1i*KX.*(data.PHI(:,:,iz,itu_)))';
+    ni(:,:,iz) = ifourier_GENE(data.DENS_I(:,:,iz,itu_))';
 end
 
 %% FILM options
@@ -92,9 +93,10 @@ while(t_<Tfin && it <= Nstep)
     if Evolve_U && (itu_old ~= itu_)
         % updating the velocity field and density field
         for iz = 1:data.Nz
-            Ux(:,:,iz) = ifourier_GENE( 1i*KY.*(data.PHI(:,:,iz,itu_)))'*sqrt(data.scale);
-            Uy(:,:,iz) = ifourier_GENE(-1i*KX.*(data.PHI(:,:,iz,itu_)))'*sqrt(data.scale);
-            ni(:,:,iz) = ifourier_GENE(data.DENS_I(:,:,iz,itu_))'*sqrt(data.scale);
+            Ux(:,:,iz) = ifourier_GENE( 1i*KY.*(data.PHI(:,:,iz,itu_)))';
+            Uy(:,:,iz) = ifourier_GENE(-1i*KX.*(data.PHI(:,:,iz,itu_)))';
+%             ni(:,:,iz) = ifourier_GENE(data.DENS_I(:,:,iz,itu_))';
+            ni(:,:,iz) = ifourier_GENE(data.Ni00(:,:,iz,itu_))';
         end
     end
     % evolve each tracer
@@ -210,10 +212,19 @@ while(t_<Tfin && it <= Nstep)
     if SHOW_FILM && (~Evolve_U || (itu_old ~= itu_))
     % updating the velocity field
         clf(fig);
-        F2P = ifourier_GENE(data.PHI(:,:,iz,itu_))';
+        switch field2plot
+            case 'phi'
+                F2P = ifourier_GENE(data.PHI(:,:,iz,itu_))';
+            case 'ni'
+                 F2P = ifourier_GENE(data.DENS_I(:,:,iz,itu_))';
+            case 'Ni00'
+                 F2P = ifourier_GENE(data.Ni00(:,:,iz,itu_))';
+            case 'none'
+                 F2P = 0*XX_;
+        end
         scale = max(max(abs(F2P))); % Scaling to normalize
         pclr = pcolor(XX_,YY_,F2P/scale); 
-        colormap(bluewhitered);
+        colormap(bluewhitered); caxis([-1 1]);
         set(pclr, 'edgecolor','none'); hold on; caxis([-2,2]); shading interp
         for ip = 1:Np
             ia0 = max(1,it-Na);
diff --git a/matlab/load/load_gene_data.m b/matlab/load/load_gene_data.m
index f3ae8237..4ef290a3 100644
--- a/matlab/load/load_gene_data.m
+++ b/matlab/load/load_gene_data.m
@@ -47,15 +47,15 @@ momfile = 'mom_ions.dat.h5';
 for jt = 1:numel(DATA.Ts3D)
     t = DATA.Ts3D(jt);
     [~, it] = min(abs(DATA.Ts3D-t));
-
-    tmp = h5read([folder,momfile],['/mom_ions/dens/',sprintf('%10.10d',it-1)]);
-    DATA.DENS_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
-
-    tmp = h5read([folder,momfile],['/mom_ions/T_par/',sprintf('%10.10d',it-1)]);
-    DATA.TPAR_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
- 
-    tmp = h5read([folder,momfile],['/mom_ions/T_perp/',sprintf('%10.10d',it-1)]);
-    DATA.TPER_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
+% 
+%     tmp = h5read([folder,momfile],['/mom_ions/dens/',sprintf('%10.10d',it-1)]);
+%     DATA.DENS_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
+% 
+%     tmp = h5read([folder,momfile],['/mom_ions/T_par/',sprintf('%10.10d',it-1)]);
+%     DATA.TPAR_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
+%  
+%     tmp = h5read([folder,momfile],['/mom_ions/T_perp/',sprintf('%10.10d',it-1)]);
+%     DATA.TPER_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
     
     tmp = h5read([folder,phifile],['/field/phi/',sprintf('%10.10d',it-1)]);
     DATA.PHI(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
diff --git a/matlab/load/quick_gene_load.m b/matlab/load/quick_gene_load.m
index 015406b7..0e0e991d 100644
--- a/matlab/load/quick_gene_load.m
+++ b/matlab/load/quick_gene_load.m
@@ -41,7 +41,7 @@
 % fname ='GENE_LIN_Kn_1.6_KT_0.4_nu_0_32x16.txt';
 % fname ='GENE_LIN_Kn_2.5_KT_0.625_nu_0_32x16.txt';
 path = '/home/ahoffman/gene/linear_CBC_results/';
-fname = 'CBC_100_ky_1e-1_to_1e1.txt';
+fname = 'CBC_linear.txt';
 data_ = load([path,fname]);
 
 figure
diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m
index b34731e2..f4deff2a 100644
--- a/matlab/plot/plot_ballooning.m
+++ b/matlab/plot/plot_ballooning.m
@@ -3,38 +3,53 @@ function [FIG] = plot_ballooning(data,options)
     [~,it0] = min(abs(data.Ts3D - options.time_2_plot(1)));
     [~,it1] = min(abs(data.Ts3D - options.time_2_plot(end)));
     [~,ikyarray] = min(abs(data.ky - options.kymodes));
-    phi_real=mean(real(data.PHI(:,:,:,it0:it1)),4);
-    phi_imag=mean(imag(data.PHI(:,:,:,it0:it1)),4);
-    
+%     phi_real=mean(real(data.PHI(:,:,:,it0:it1)),4);
+%     phi_imag=mean(imag(data.PHI(:,:,:,it0:it1)),4);
+    phi_real=real(data.PHI(:,:,:,it1));
+    phi_imag=imag(data.PHI(:,:,:,it1));
     % Apply baollooning tranform
     for iky=ikyarray
         dims = size(phi_real);
-
-        if data.SHEAR > 0
-            idx=[0:data.Nkx/2 -(data.Nkx/2-1):-1];
-            ikxlim = dims(2);
+        Nkx  = dims(2);
+        is   = max(1,iky-1);
+        Npi  = (Nkx-1)-2*(is-1);
+        if(Npi <= 0)
+            break
+        elseif(Npi == 1)
+            ordered_ikx = 1;
         else
-            idx = 0;
-            ikxlim = 1;
+            tmp_ = (Nkx-is+1):-is:(Nkx/2+2);
+            ordered_ikx = [tmp_(end:-1:1), 1:is:(Nkx/2)];
         end
+
+        idx=data.kx./data.kx(2);
+        idx= idx(ordered_ikx);
+        Nkx = numel(idx);
+
+        phib_real = zeros(  Nkx*dims(3)  ,1);
+        phib_imag = phib_real;
+        b_angle   = phib_real;
         
-        phib_real = zeros(  ikxlim*dims(3)  ,1);
-        phib_imag= phib_real;
-        b_angle = phib_real;
-        
-        for ikx =1:ikxlim
-            start_ =  (ikx -1)*dims(3) +1;
-            end_ =  ikx*dims(3);
+        kk_=[];
+        ss_=[];
+        ee_=[];
+        for i_ =1:Nkx
+            start_ =  (i_ -1)*dims(3) +1;
+            end_ =  i_*dims(3);
+            ikx  = ordered_ikx(i_);
             phib_real(start_:end_)  = phi_real(iky,ikx,:);
             phib_imag(start_:end_)  = phi_imag(iky,ikx,:);
+            kk_ = [kk_ data.kx(ikx)];
+            ss_ = [ss_ start_];
+            ee_ = [ee_ end_];
         end
 
         % Define ballooning angle
-        Nkx = numel(data.kx)-1; coordz = data.z;
-        for ikx =1: ikxlim
+        coordz = data.z;
+        for i_ =1: Nkx
             for iz=1:dims(3)
-                ii = dims(3)*(ikx -1) + iz;
-                b_angle(ii) =coordz(iz) + 2*pi*idx(ikx);
+                ii = dims(3)*(i_ -1) + iz;
+                b_angle(ii) =coordz(iz) + 2*pi*idx(i_)./is;
             end
         end
         
diff --git a/matlab/plot/spectrum_1D.m b/matlab/plot/spectrum_1D.m
index 8c9c90b2..14ebd494 100644
--- a/matlab/plot/spectrum_1D.m
+++ b/matlab/plot/spectrum_1D.m
@@ -31,6 +31,9 @@ switch options.NAME
     case 'n_e'
     FIGURE.fig = figure; FIGURE.FIGNAME = ['ne_spectrum_',data.PARAMS]; 
     fieldname = 'elec. dens.';
+    case 'N_i^{00}'
+    FIGURE.fig = figure; FIGURE.FIGNAME = ['Ni00_spectrum_',data.PARAMS]; 
+    fieldname = 'ion gdens.';
 end
 
 PLOT2D = 0;
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index c80985ea..99ceca4d 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -46,17 +46,17 @@ options.INTERP    = 1;
 options.POLARPLOT = 0;
 % options.NAME      = '\phi';
 % options.NAME      = 'N_i^{00}';
-% options.NAME      = 'v_x';
+options.NAME      = 'v_y';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
-options.NAME      = 'n_i';
+% 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(20:2:end);
-options.TIME      = [700:1:750];
+options.TIME      =  data.Ts3D(1:20:end);
+% options.TIME      = [750:1:1000];
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 create_film(data,options,'.gif')
@@ -70,7 +70,7 @@ options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
 % options.NAME      = '\phi';
 % options.NAME      = 'n_e';
-% options.NAME      = 'N_i^{00}';
+options.NAME      = 'N_i^{00}';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
@@ -78,7 +78,7 @@ options.PLAN      = 'xy';
 % options.NAME      'f_i';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = [200 800 1500];
+options.TIME      = [500 800 1000];
 data.a = data.EPS * 2e3;
 fig = photomaton(data,options);
 % save_figure(data,fig)
@@ -131,7 +131,7 @@ end
 
 if 0
 %% Time averaged spectrum
-options.TIME   = [500 700];
+options.TIME   = [300 600];
 options.NORM   =1;
 % options.NAME   = '\phi';
 % options.NAME      = 'N_i^{00}';
@@ -139,7 +139,7 @@ options.NAME   ='\Gamma_x';
 options.PLAN   = 'kxky';
 options.COMPZ  = 'avg';
 options.OK     = 0;
-options.COMPXY = 'avg'; % avg/sum/max/zero/ 2D plot otherwise
+options.COMPXY = '2D'; % avg/sum/max/zero/ 2D plot otherwise
 options.COMPT  = 'avg';
 options.PLOT   = 'semilogy';
 fig = spectrum_1D(data,options);
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index bd2e6202..822daf31 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -11,7 +11,8 @@
 % folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
 % folder = '/misc/gene_results/LD_zpinch_1.6/';
 % folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
-folder = '/misc/gene_results/ZP_kn_2.5_large_box/';
+% folder = '/misc/gene_results/ZP_kn_2.5_large_box/';
+folder = '/misc/gene_results/CBC/128x64x16x24x12/';
 gene_data = load_gene_data(folder);
 gene_data = invert_kxky_to_kykx_gene_results(gene_data);
 if 1
@@ -23,7 +24,7 @@ options.ST_FIELD = '\phi';          % chose your field to plot in spacetime diag
 options.INTERP   = 1;
 gene_data.FIGDIR = folder;
 fig = plot_radial_transport_and_spacetime(gene_data,options);
-save_figure(gene_data,fig)
+save_figure(gene_data,fig,'.png')
 end
 
 if 0
@@ -39,8 +40,8 @@ options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
 % options.NAME      = 'Q_x';
-% options.NAME      = '\phi';
-options.NAME      = 'n_i';
+options.NAME      = '\phi';
+% options.NAME      = 'n_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
 options.PLAN      = 'xy';
@@ -56,13 +57,13 @@ end
 if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Options
-options.INTERP    = 1;
+options.INTERP    = 1;                  
 options.POLARPLOT = 0;
-% options.NAME      = '\phi';
+options.NAME      = '\phi';
 % options.NAME      = 'v_y';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
-options.NAME      = 'n_i';
+% options.NAME      = 'n_i';
 options.PLAN      = 'kxky';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
diff --git a/wk/g2k3_transport_scaling.m b/wk/g2k3_transport_scaling.m
index 3f06c1a2..c75cd198 100644
--- a/wk/g2k3_transport_scaling.m
+++ b/wk/g2k3_transport_scaling.m
@@ -29,10 +29,10 @@ figure;
 % plot(K_N,G_LR); hold on;
 % plot(K_N,G_LD); hold on;
 % plot(K_N,G_CL,'--k'); hold on;
-k = 0.35;
-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;
+k = 0.3;
+plot(K_N,g_SG.^2/k.^3./K_N); hold on;
+plot(K_N,g_DG.^2/k.^3./K_N); hold on;
+plot(K_N,g_LR.^2/k.^3./K_N); hold on;
+plot(K_N,g_LD.^2/k.^3./K_N); 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 ef37d89e..7088014c 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -28,11 +28,12 @@ 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_256x128x5x3_new_NL';
-outfile ='Zpinch_rerun/Kn_2.5_256x64x5x3';
+% outfile ='Zpinch_rerun/Kn_2.5_256x128x5x3';
+% outfile ='Zpinch_rerun/Kn_2.5_312x196x5x3_Lx_400_Ly_200';
+% outfile ='Zpinch_rerun/Kn_2.5_256x64x5x3';
 % outfile ='Zpinch_rerun/Kn_2.0_200x48x9x5_large_box';
 % outfile ='Zpinch_rerun/Kn_2.0_256x64x9x5_Lx_240_Ly_120';
 % outfile ='Zpinch_rerun/Kn_1.6_256x128x7x4';
 % outfile ='Zpinch_rerun/Kn_1.6_200x48x11x6';
-% outfile ='Zpinch_rerun/Kn_1.6_200x48x21x11';
+outfile ='Zpinch_rerun/Kn_1.6_256x128x21x11';
 run analysis_HeLaZ
diff --git a/wk/quick_run.m b/wk/quick_run.m
index 63c31a65..737428b0 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -13,23 +13,23 @@ EXECNAME = 'helaz3';
 CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 0.5;   % Collision frequency
+NU      = 0.01;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
-K_N     = 6.0;%2.0;   % Density gradient drive
-K_T     = 0;%0.25*K_N;   % Temperature '''
+K_N     = 2.22;%2.0;   % Density gradient drive
+K_T     = 6.92;%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)
-KIN_E   = 1;     % 1: kinetic electrons, 2: adiabatic electrons
+KIN_E   = 0;     % 1: kinetic electrons, 2: adiabatic electrons
 %% GRID PARAMETERS
 PMAXE   = 4;     % Hermite basis size of electrons
 JMAXE   = 2;     % Laguerre "
 PMAXI   = 4;     % " ions
 JMAXI   = 2;     % "
-NX      = 32;    % real space x-gridpoints
-NY      = 1;     %     ''     y-gridpoints
-LX      = 100;   % Size of the squared frequency domain
-LY      = 2*pi/0.5;     % Size of the squared frequency domain
+NX      = 12;    % real space x-gridpoints
+NY      = 8;     %     ''     y-gridpoints
+LX      = 2*pi/0.1;   % Size of the squared frequency domain
+LY      = 2*pi/0.1;     % Size of the squared frequency domain
 NZ      = 16;     % number of perpendicular planes (parallel grid)
 NPOL    = 1;
 SG      = 0;     % Staggered z grids option
@@ -37,11 +37,11 @@ SG      = 0;     % Staggered z grids option
 % GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
 GEOMETRY= 's-alpha';
 Q0      = 1.4;    % safety factor
-SHEAR   = 0.796;    % magnetic shear (Not implemented yet)
+SHEAR   = 0.8;    % magnetic shear (Not implemented yet)
 EPS     = 0.18;    % inverse aspect ratio
 %% TIME PARMETERS
-TMAX    = 50;  % Maximal time unit
-DT      = 1e-3;   % Time step
+TMAX    = 20;  % Maximal time unit
+DT      = 5e-3;   % Time step
 SPS0D   = 1;      % Sampling per time unit for 2D arrays
 SPS2D   = 0;      % Sampling per time unit for 2D arrays
 SPS3D   = 1;      % Sampling per time unit for 2D arrays
@@ -49,11 +49,11 @@ SPS5D   = 1/5;    % Sampling per time unit for 5D arrays
 SPSCP   = 0;    % Sampling per time unit for checkpoints
 JOB2LOAD= -1;
 %% OPTIONS
-SIMID   = 'shear_testcase_Pan22_linear';  % Name of the simulation
+SIMID   = 'linear_CBC';  % 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      = 'SG';
+CO      = 'DG';
 GKCO    = 0; % gyrokinetic operator
 ABCO    = 1; % interspecies collisions
 INIT_ZF = 0; ZF_AMP = 0.0;
@@ -77,23 +77,24 @@ INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
 MU_X    = MU;     %
 MU_Y    = MU;     %
 N_HD    = 4;
-MU_Z    = 0.05;     %
+MU_Z    = 2.0;     %
 MU_P    = 0.0;     %
 MU_J    = 0.0;     %
 LAMBDAD = 0.0;
-NOISE0  = 1.0e-5; % Init noise amplitude
-BCKGD0  = 0.0;    % Init background
+NOISE0  = 0.0e-5; % Init noise amplitude
+BCKGD0  = 1.0;    % Init background
 GRADB   = 1.0;
 CURVB   = 1.0;
 %%-------------------------------------------------------------------------
 %% RUN
 setup
-system(['rm fort*.90']);
+% system(['rm fort*.90']);
 % Run linear simulation
 if RUN
 %     system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
+%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
 %     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
-    system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',HELAZDIR,'bin/',EXECNAME,' 1 1 2 0; cd ../../../wk'])
+    system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 2 1 3 0; cd ../../../wk'])
 %     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
 end
 
@@ -106,7 +107,7 @@ JOBNUMMIN = 00; JOBNUMMAX = 00;
 data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
 
 %% Short analysis
-if 1
+if 0
 %% linear growth rate (adapted for 2D zpinch and fluxtube)
 trange = [0.5 1]*data.Ts3D(end);
 nplots = 1;
@@ -117,10 +118,10 @@ 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
+if 1
 %% Ballooning plot
-options.time_2_plot = [0.5 0.5]*data.Ts3D(end);
-options.kymodes     = [0.5];
+options.time_2_plot = [1]*data.Ts3D(end);
+options.kymodes     = [0.1 0.2 0.3];
 options.normalized  = 1;
 options.field       = 'phi';
 fig = plot_ballooning(data,options);
@@ -159,9 +160,9 @@ if 0
 %% Mode evolution
 options.NORMALIZED = 1;
 options.K2PLOT = 1;
-options.TIME   = [0.8 1]*data.Ts3D(end);
+options.TIME   = [0 1]*data.Ts3D(end);
 options.NMA    = 1;
-options.NMODES = 2;
+options.NMODES = 1;
 options.iz     = 9;
 fig = mode_growth_meter(data,options);
 save_figure(gbms_dat,fig)
-- 
GitLab