scripts save

matlab/evolve_tracers.m

 % Options
-SHOW_FILM = 0;
+SHOW_FILM = 1;
 field2plot  ='phi';
 INIT     = 'lin';   % lin (for a line)/ round (for a small round)/ gauss for random
 U_TIME   = 1000;     % >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     = 100;
+Tfin     = 500;
 dt_      = 0.1;
 Nstep    = ceil(Tfin/dt_);
 % Init tracers
@@ -181,8 +181,7 @@ while(t_<Tfin && it <= Nstep)
             u___  =  linterp(u__0,u__1,a__i);
 
 %             push the particle
-%             q = sign(-u___(3));
-            q = 1;%-u___(3);
+            q = sign(-u___(3));
 %             q =1;
             x_ = x_ + dt_*u___(1)*q;
             y_ = y_ + dt_*u___(2)*q;
@@ -224,19 +223,20 @@ while(t_<Tfin && it <= Nstep)
         end
         scale = max(max(abs(F2P))); % Scaling to normalize
         pclr = pcolor(XX_,YY_,F2P/scale); 
-        caxis([-1 1]); colormap(bluewhitered); colorbar;
+        caxis([-1 1]); colormap(bluewhitered); %colorbar;
         set(pclr, 'edgecolor','none'); hold on; caxis([-1+dimmed,1+dimmed]); shading interp
         for ip = 1:Np
             ia0 = max(1,it-Na);
             plot(Traj_x(ip,ia0:it),Traj_y(ip,ia0:it),'.','Color',color(ip,:)); hold on
         end
         for ip = 1:Np
-            plot(Traj_x(ip, 1),Traj_y(ip, 1),'xk'); hold on
+%             plot(Traj_x(ip, 1),Traj_y(ip, 1),'xk'); hold on
             plot(Traj_x(ip,it),Traj_y(ip,it),'ok','MarkerFaceColor',color(ip,:)); hold on
         end
         title(['$t \approx$', sprintf('%.3d',ceil(data.Ts3D(itu_)))]);
         axis equal
         xlim([xmin xmax]); ylim([ymin ymax]);
+        xlabel('$x/\rho_s$'); ylabel('$y/\rho_s$');
         drawnow
         % Capture the plot as an image 
         frame = getframe(fig); 
@@ -303,10 +303,14 @@ xlim(U_TIME + [0 Tfin]);
 subplot(222);
     itf = floor(Nt/2); %fit end time
     % x^2 displacement
-    plot(time_,mean(xtot.^2,1),'DisplayName','$\langle x.^2\rangle_p$'); hold on
+    x2tot = mean(xtot.^2,1);
+    plot(time_-time_(1),x2tot,'DisplayName','$\langle x.^2\rangle_p$'); hold on
     fit = polyfit(time_(1:itf),mean(xtot(:,1:itf).^2,1),1);
-    plot(time_,fit(1)*time_+fit(2),'--k'); hold on
-    ylabel('$\langle x^2 \rangle_p$');
+    plot(time_-time_(1),fit(1)*time_+fit(2),'--k'); hold on
+    % Put a linear growth from the first point to check if it super/sub
+    % diffusive
+%     loglog(time_-time_(1),(time_-time_(1))*x2tot(2)/(time_(2)-time_(1)),'--k'); hold on
+%     ylabel('$\langle x^2 \rangle_p$');
 
 %     % y^2 displacement
 %     fit = polyfit(time_(1:itf),mean(ytot(:,1:itf).^2,1),1);

matlab/extract_fig_data.m

@@ -4,7 +4,7 @@
 % tw = [3000 4000];
 % tw = [4000 4500];
 % tw = [4500 5000];
-tw = [500 1000];
+tw = [00 5000];
 
 fig = gcf;
 axObjs = fig.Children;
@@ -19,13 +19,13 @@ end
 % n1 = numel(X_);
 [~,n0] = min(abs(X_-tw(1)));
 [~,n1] = min(abs(X_-tw(2)));
-mvm = @(x) movmean(x,50);
-shift = X_(n0);
+mvm = @(x) movmean(x,1);
+shift = 0;%X_(n0);
 % shift = 0;
+skip = 50;
 
 figure;
-plot(X_(n0:end),Y_(n0:end));
-plot(mvm(X_(n0:n1)-shift),mvm(Y_(n0:n1))); hold on;
+plot(mvm(X_(n0:skip:n1)-shift),mvm(Y_(n0:skip:n1))); hold on;
 
 % t0 = ceil(numel(X_)*0.2); t1 = numel(X_);
 avg= mean(Y_(n0:n1)); dev = std(Y_(n0:n1));

wk/analysis_gene.m

@@ -12,9 +12,9 @@ addpath(genpath([helazdir,'matlab/load'])) % ... add
 % 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/Z-pinch/HP_fig_2a_mu_1e-2/';
+% folder = '/misc/gene_results/Z-pinch/HP_fig_2b_mu_5e-2/';
+% folder = '/misc/gene_results/Z-pinch/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_HP_kn_1.6_nuv_3.2/';
@@ -29,7 +29,7 @@ addpath(genpath([helazdir,'matlab/load'])) % ... add
 % folder = '/misc/gene_results/CBC/KT_13_large_box_128x64x16x24x12/';
 % folder = '/misc/gene_results/CBC/Lapillone_Fig6/';
 % folder = '/misc/gene_results/Z-pinch/HP_kN_1.6_adapt_mu_01/';
-folder = '/misc/gene_results/miller/';
+% folder = '/misc/gene_results/miller/';
 gene_data = load_gene_data(folder);
 gene_data = invert_kxky_to_kykx_gene_results(gene_data);
 if 1
@@ -76,9 +76,9 @@ if 0
 % Options
 options.INTERP    = 1;
 options.POLARPLOT = 0;
-options.NAME      = '\phi';
+% options.NAME      = '\phi';
 % options.NAME      = 'v_y';
-% options.NAME      = 'n_i^{NZ}';
+options.NAME      = 'n_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
 options.PLAN      = 'xy';
@@ -131,15 +131,15 @@ end
 
 if 0
 %% Time averaged spectrum
-options.TIME   = 300:600;
+options.TIME   = [4000 8000];
 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;
-options.COMPXY = 'zero';
+options.COMPXY = 'avg'; % avg/sum/max/zero/ 2D plot otherwise
 options.COMPT  = 'avg';
 options.PLOT   = 'semilogy';
 fig = spectrum_1D(gene_data,options);

wk/analysis_gyacomo.m

@@ -34,7 +34,7 @@ disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))])
 options.TAVG_0   = data.TJOB_SE(i_);%0.4*data.Ts3D(end);
 options.TAVG_1   = data.TJOB_SE(i_+1);%0.9*data.Ts3D(end); % Averaging times duration
 options.NCUT     = 4;              % Number of cuts for averaging and error estimation
-options.NMVA     = 1;              % Moving average for time traces
+options.NMVA     = 100;              % 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 = '\phi';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
 options.INTERP   = 0;
@@ -59,13 +59,13 @@ if 0
 % Options
 options.INTERP    = 1;
 options.POLARPLOT = 0;
-options.NAME      = '\phi';
+% options.NAME      = '\phi';
 % options.NAME      = '\omega_z';
 % options.NAME      = 'N_i^{00}';
 % 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';
@@ -154,9 +154,9 @@ if 0
 %% Time averaged spectrum
 options.TIME   = [2000 3000];
 options.NORM   =1;
-options.NAME   = '\phi';
+% options.NAME   = '\phi';
 % options.NAME      = 'N_i^{00}';
-% options.NAME   ='\Gamma_x';
+options.NAME   ='\Gamma_x';
 options.PLAN   = 'kxky';
 options.COMPZ  = 'avg';
 options.OK     = 0;

wk/header_2DZP_results.m

@@ -200,10 +200,25 @@ resdir ='';
 %% Convergence cases kN = (1.6 2.2) nu = (0.01 1.0)
 % resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x3x2_mu_0.01';
 % resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x5x3_mu_0.01';
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x7x4_mu_0.01';
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x9x5_mu_0.01';
 
 % resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x3x2_mu_0.01';
-resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x7x4_mu_0.01';
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x5x3_mu_0.01';
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x7x4_mu_0.01';
+resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x9x5_mu_0.01';
 
-JOBNUMMIN = 00; JOBNUMMAX = 10;
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x3x2_nu_0.1';
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x5x3_nu_0.1';
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x7x4_nu_0.1';
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x9x5_nu_0.1';
+
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x3x2_nu_0.1';
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x5x3_nu_0.1';
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x7x4_nu_0.1';
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x9x5_nu_0.1';
+
+
+JOBNUMMIN = 00; JOBNUMMAX = 03;
 resdir = ['results/',resdir];
 run analysis_gyacomo
\ No newline at end of file

wk/lin_EPY.m

@@ -4,21 +4,21 @@
 % for benchmark and debugging purpose since it makes matlab "busy"
 %
 SIMID   = 'lin_EPY';  % Name of the simulation
-SIMID   = 'dbg';  % Name of the simulation
-RUN     = 1; % To run or just to load
+% SIMID   = 'dbg';  % Name of the simulation
+RUN     = 0; % To run or just to load
 addpath(genpath('../matlab')) % ... add
 default_plots_options
 PROGDIR  = '/home/ahoffman/gyacomo/';
-EXECNAME = 'gyacomo_dbg';
-% EXECNAME = 'gyacomo';
+% EXECNAME = 'gyacomo_dbg';
+EXECNAME = 'gyacomo';
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Set Up parameters
 CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 0.01;           % Collision frequency
+NU      = 0.1;           % Collision frequency
 TAU     = 1.0;            % e/i temperature ratio
-K_Ne    = 2.2;            % ele Density '''
+K_Ne    = 1.6;            % ele Density '''
 K_Te    = K_Ne/4;            % ele Temperature '''
 K_Ni    = K_Ne;            % ion Density gradient drive
 K_Ti    = K_Ni/4;            % ion Temperature '''
@@ -33,7 +33,7 @@ JMAXE   = J;     % Laguerre "
 PMAXI   = P;     % " ions
 JMAXI   = J;     % "
 NX      = 2;    % real space x-gridpoints
-NY      = 2;     %     ''     y-gridpoints
+NY      = 100;     %     ''     y-gridpoints
 LX      = 2*pi/0.8;   % Size of the squared frequency domain
 LY      = 120;%2*pi/0.05;     % Size of the squared frequency domain
 NZ      = 1;    % number of perpendicular planes (parallel grid)
@@ -46,11 +46,11 @@ SHEAR   = 0.8;    % magnetic shear
 NEXC    = 1;      % To extend Lx if needed (Lx = Nexc/(kymin*shear))
 EPS     = 0.18;   % inverse aspect ratio
 %% TIME PARMETERS
-TMAX    = 200;  % Maximal time unit
+TMAX    = 500;  % Maximal time unit
 DT      = 1e-2;   % Time step
 SPS0D   = 1;      % Sampling per time unit for 2D arrays
 SPS2D   = 0;      % Sampling per time unit for 2D arrays
-SPS3D   = 2;      % Sampling per time unit for 2D arrays
+SPS3D   = 1/5;      % Sampling per time unit for 2D arrays
 SPS5D   = 1/5;    % Sampling per time unit for 5D arrays
 SPSCP   = 0;    % Sampling per time unit for checkpoints
 JOB2LOAD= -1;
@@ -111,7 +111,7 @@ data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from
 if 1
 %% linear growth rate (adapted for 2D zpinch and fluxtube)
 options.TRANGE = [0.5 1]*data.Ts3D(end);
-options.NPLOTS = 1; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
+options.NPLOTS = 2; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
 options.GOK    = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
 lg = compute_fluxtube_growth_rate(data,options);
 [gmax,     kmax] = max(lg.g_ky(:,end));
@@ -157,13 +157,13 @@ options.kzky = 0;
 [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
 save_figure(data,fig)
 end
-if 0
+if 1
 %% Mode evolution
 options.NORMALIZED = 0;
-options.K2PLOT = 1;
+options.K2PLOT = 10;
 options.TIME   = [0:1000];
 options.NMA    = 1;
-options.NMODES = 1;
+options.NMODES = 10;
 options.iz     = 'avg';
 fig = mode_growth_meter(data,options);
 save_figure(data,fig,'.png')