diff --git a/matlab/evolve_tracers.m b/matlab/evolve_tracers.m
index 34df7362df15d32c724bb8be77521a9d3eeb3445..4725ede3741e371fce409a6c9c85fb5398ec46c5 100644
--- a/matlab/evolve_tracers.m
+++ b/matlab/evolve_tracers.m
@@ -1,10 +1,10 @@
% 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);
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index b9d69af104a76dfbf926fd4f59ebef3161a25934..a036d18452ace4fd9aaca42f567028ae996b8663 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 = [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));
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index de09d1690e7bb9c26fee4bafd7b35ed833e0cc4e..53f11c4f9f6cb9dffb0aa4291a7d347b1c5b0f66 100644
--- a/wk/analysis_gene.m
+++ b/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);
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index d521f606b1347651fa70ee0a74f87382961ab564..2a875f63bc57ccab9a3aacc780124248ef6bcd97 100644
--- a/wk/analysis_gyacomo.m
+++ b/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;
diff --git a/wk/header_2DZP_results.m b/wk/header_2DZP_results.m
index 32653c85cff97c84df96cf0d0b3065e3e29a2d8d..cdd4074d3972b6cd85e0275ae582dea29129c9dc 100644
--- a/wk/header_2DZP_results.m
+++ b/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
diff --git a/wk/lin_EPY.m b/wk/lin_EPY.m
index b49201d025f2e9f4872fcc826be6d1d0040bc6c0..2ddd756a47e22d4e3b808e4f0c8dd019d6da95fd 100644
--- a/wk/lin_EPY.m
+++ b/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')