diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 35c794914a7af795f2df5f97827d27e97872fb9e..3ee858a779b539957b602fd00cbb197fdb1a4f68 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -1,4 +1,4 @@
-function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS,CODE)
+function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
%Compute steady radial transport
tend = OPTIONS.TAVG_1; tstart = OPTIONS.TAVG_0;
[~,its0D] = min(abs(DATA.Ts0D-tstart));
@@ -11,11 +11,11 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS,CODE)
Qx_infty_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE;
Qx_infty_std = std (DATA.HFLUX_X(its0D:ite0D))*SCALE;
% disp(['Q_x=',sprintf('%2.2e',Qx_infty_avg),'+-',sprintf('%2.2e',Qx_infty_std)]);
- f_avg_z = squeeze(mean(DATA.PHI(:,:,:,:),3));
- [~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3)));
- ikzf = min([ikzf,DATA.Nky]);
- Ns3D = numel(DATA.Ts3D);
- [KX, KY] = meshgrid(DATA.kx, DATA.ky);
+% f_avg_z = squeeze(mean(DATA.PHI(:,:,:,:),3));
+% [~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3)));
+% ikzf = min([ikzf,DATA.Nky]);
+% Ns3D = numel(DATA.Ts3D);
+% [KX, KY] = meshgrid(DATA.kx, DATA.ky);
%% error estimation
DT_ = (tend-tstart)/OPTIONS.NCUT;
Qx_ee = zeros(1,OPTIONS.NCUT);
@@ -27,36 +27,36 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS,CODE)
Qx_avg = mean(Qx_ee);
Qx_err = std(Qx_ee);
disp(['Q_avg=',sprintf('%2.2e',Qx_avg),'+-',sprintf('%2.2e',Qx_err)]);
- %% computations
-
- % Compute zonal and non zonal energies
- 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(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
- % Compute thermodynamic entropy Eq.(5) Navarro et al. 2012 PoP
- % 1/2 sum_p sum_j Napj^2(k=0) (avg z)
- Nipjz = sum(sum(sum(sum(conj(DATA.Nipj).*DATA.Nipj))));
- ff = trapz(DATA.z,Nipjz,5);
- E_TE = 0.5*squeeze(ff);
- % Compute electrostatic energy
- E_ES = zeros(size(DATA.Ts5D));
- bi = sqrt(KX.^2+KY.^2)*DATA.sigma(1)*sqrt(2*DATA.tau(1)); %argument of the kernel
- for it5D = 1:numel(DATA.Ts5D)
- [~,it3D] = min(abs(DATA.Ts3D-DATA.Ts5D(it5D)));
- for in = 1:DATA.Jmaxi
- Knphi = kernel(in-1,bi).*squeeze(trapz(DATA.z,DATA.PHI(:,:,:,it3D),3));
- Ni0n_z= squeeze(trapz(DATA.z,DATA.Nipj(1,in,:,:,:,it5D),5));
- E_ES(it5D) = 0.5*sum(sum(abs(conj(Knphi).*Ni0n_z)));
- end
- end
+% %% computations
+%
+% % Compute zonal and non zonal energies
+% 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(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
+% % Compute thermodynamic entropy Eq.(5) Navarro et al. 2012 PoP
+% % 1/2 sum_p sum_j Napj^2(k=0) (avg z)
+% Nipjz = sum(sum(sum(sum(conj(DATA.Nipj).*DATA.Nipj))));
+% ff = trapz(DATA.z,Nipjz,5);
+% E_TE = 0.5*squeeze(ff);
+% % Compute electrostatic energy
+% E_ES = zeros(size(DATA.Ts5D));
+% bi = sqrt(KX.^2+KY.^2)*DATA.sigma(1)*sqrt(2*DATA.tau(1)); %argument of the kernel
+% for it5D = 1:numel(DATA.Ts5D)
+% [~,it3D] = min(abs(DATA.Ts3D-DATA.Ts5D(it5D)));
+% for in = 1:DATA.Jmaxi
+% Knphi = kernel(in-1,bi).*squeeze(trapz(DATA.z,DATA.PHI(:,:,:,it3D),3));
+% Ni0n_z= squeeze(trapz(DATA.z,DATA.Nipj(1,in,:,:,:,it5D),5));
+% E_ES(it5D) = 0.5*sum(sum(abs(conj(Knphi).*Ni0n_z)));
+% end
+% end
%% Figure
clr_ = lines(20);
mvm = @(x) movmean(x,OPTIONS.NMVA);
FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; %set(gcf, 'Position', [500, 1000, 1000, 600])
- FIGURE.ax1 = subplot(3,1,1,'parent',FIGURE.fig);
+ FIGURE.ax1 = subplot(2,1,1,'parent',FIGURE.fig);
plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'--',...
'color',clr_((DATA.Pmaxi-1)/2,:),...
'DisplayName',['$\Gamma_x$ ',DATA.paramshort]); hold on;
@@ -80,16 +80,16 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
title({DATA.param_title,...
['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$, Q^{\infty} = $',num2str(Qx_infty_avg)]});
- %% Free energy
- FIGURE.ax2 = subplot(3,1,2,'parent',FIGURE.fig);
- yyaxis left
- plot(DATA.Ts5D,E_TE,'DisplayName','$\epsilon_f$'); hold on;
- ylabel('Entropy');%('$\epsilon_f$')
- yyaxis right
- plot(DATA.Ts5D,E_ES,'DisplayName','$\epsilon_\phi$');
- ylabel('ES energy');%('$\epsilon_\phi$')
- xlim([DATA.Ts5D(1), DATA.Ts5D(end)]);
- xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
+% %% Free energy
+% FIGURE.ax2 = subplot(3,1,2,'parent',FIGURE.fig);
+% yyaxis left
+% plot(DATA.Ts5D,E_TE,'DisplayName','$\epsilon_f$'); hold on;
+% ylabel('Entropy');%('$\epsilon_f$')
+% yyaxis right
+% plot(DATA.Ts5D,E_ES,'DisplayName','$\epsilon_\phi$');
+% ylabel('ES energy');%('$\epsilon_\phi$')
+% xlim([DATA.Ts5D(1), DATA.Ts5D(end)]);
+% xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
%% radial shear radial profile
% computation
Ns3D = numel(DATA.Ts3D);
@@ -108,7 +108,7 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
f2plot = toplot.FIELD;
dframe = ite3D - its3D;
clim = max(max(max(abs(plt(f2plot(:,:,:))))));
- FIGURE.ax3 = subplot(3,1,3,'parent',FIGURE.fig);
+ FIGURE.ax3 = subplot(2,1,2,'parent',FIGURE.fig);
[TY,TX] = meshgrid(DATA.x,DATA.Ts3D(toplot.FRAMES));
pclr = pcolor(TX,TY,squeeze(plt(f2plot))');
set(pclr, 'edgecolor','none');
diff --git a/matlab/profiler.m b/matlab/profiler.m
index d509975b60c9b4919c8e5cb11c3c5d00acaab895..5056f61e4e709e2f65f094e9dd41484c57c8bfb0 100644
--- a/matlab/profiler.m
+++ b/matlab/profiler.m
@@ -95,7 +95,7 @@ end
p1 = area(xx_,yy_,'LineStyle','none');
for i = 1:N_T; p1(i).FaceColor = colors(i,:);
% LEGEND{i} = sprintf('%s t=%1.1e[s] (%0.1f %s)',names{i},Ts_A(i),Ts_A(i)/total_Tc(end)*100,'\%');
- LEGEND{i} = [names{i},' $\hat t=$',sprintf('%1.1e[s] (%0.1f %s)',Ts_A(i),Ts_A(i)/total_Ta*100,'\%')];
+ LEGEND{i} = [names{i},' $\hat t=$',sprintf('%1.1e[s] (%0.1f %s)',Ts_A(i)/NSTEP_PER_SAMP,Ts_A(i)/total_Ta*100,'\%')];
end;
legend(LEGEND,'Location','bestoutside');
% legend('Compute RHS','Adv. fields','ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Process', 'Missing')
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index c098bc27d2332000619cd0e2cf2128fdd14720f4..cb18516afc17a82032d8d52fbd0e41a64b4d6a74 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -80,7 +80,7 @@ options.NAME = '\phi';
% options.NAME = 'Q_x';
% options.NAME = 'n_i';
% options.NAME = 'n_i-n_e';
-options.PLAN = 'xy';
+options.PLAN = 'kxky';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
@@ -105,16 +105,15 @@ options.NORMALIZE = 0;
% options.NAME = '\omega_z';
% options.NAME = 'T_i';
% options.NAME = 'n_i';
-% options.NAME = '\phi^{NZ}';
-options.NAME = 'N_i^{00}';
+options.NAME = '\phi^{NZ}';
+% options.NAME = 'N_i^{00}';
% options.NAME = 'N_i^{00}-N_e^{00}';
% options.NAME = 's_{Ex}';
% options.NAME = 'Q_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'kxky';
+options.PLAN = 'xy';
options.COMP = 'avg';
-% options.TIME = [49 50 51];
-options.TIME = data.Ts3D(51:55);
+options.TIME = [80 200];
options.RESOLUTION = 256;
data.a = data.EPS * 2e3;
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index 83f4ff9e078d4f58cae1062549b196d44342f3ce..fdd59170dabe410d0bf25d2688fbbb321e5deb8b 100644
--- a/wk/fast_analysis.m
+++ b/wk/fast_analysis.m
@@ -4,17 +4,50 @@ PARTITION = '/misc/gyacomo23_outputs/';
% PARTITION = gyacomodir;
%% CBC
-resdir = 'paper_2_GYAC23/CBC/7x4x192x96x32';
+% resdir = 'paper_2_GYAC23/CBC/7x4x192x96x32_nu_0.05_muxy_0.5_muz_0.2';
+resdir = 'paper_2_GYAC23/CBC/7x4x192x96x32_nu_0.05_muxy_1.0_muz_1.0;'
%%
-JOBNUMMIN = 00; JOBNUMMAX = 10;
+J0 = 00; J1 = 10;
-%%
+% Load basic info (grids and time traces)
DATADIR = [PARTITION,resdir,'/'];
data = {};
-data = compile_results_low_mem(data,DATADIR,JOBNUMMIN,JOBNUMMAX);
+data = compile_results_low_mem(data,DATADIR,J0,J1);
+%% Plot transport and phi radial profile
+[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
-%%
-figure;
-plot(data.Ts0D,data.HFLUX_X);
\ No newline at end of file
+options.TAVG_0 = 100;
+options.TAVG_1 = 1000;
+options.NCUT = 5; % Number of cuts for averaging and error estimation
+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 = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+options.INTERP = 0;
+options.RESOLUTION = 256;
+fig = plot_radial_transport_and_spacetime(data,options);
+
+%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Options
+options.INTERP = 0;
+options.POLARPLOT = 0;
+options.NAME = '\phi';
+% options.NAME = '\omega_z';
+% options.NAME = 'N_i^{00}';
+% options.NAME = 's_{Ey}';
+% options.NAME = 'n_i^{NZ}';
+% options.NAME = 'Q_x';
+% options.NAME = 'n_i';
+% options.NAME = 'n_i-n_e';
+options.PLAN = 'xz';
+% options.NAME = 'f_i';
+% options.PLAN = 'sx';
+options.COMP = 'avg';
+% options.TIME = data.Ts5D(end-30:end);
+% options.TIME = data.Ts3D;
+options.TIME = [0:1500];
+data.EPS = 0.1;
+data.a = data.EPS * 2000;
+options.RESOLUTION = 256;
+create_film(data,options,'.gif')