Skip to content
Snippets Groups Projects
Commit 11ff59c4 authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

ud

parent 351e5947
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
matlab/plot/plot_radial_transport_and_spacetime.m
+ 43
− 43
View file @ 11ff59c4
function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS,CODE) function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
%Compute steady radial transport %Compute steady radial transport
tend = OPTIONS.TAVG_1; tstart = OPTIONS.TAVG_0; tend = OPTIONS.TAVG_1; tstart = OPTIONS.TAVG_0;
[~,its0D] = min(abs(DATA.Ts0D-tstart)); [~,its0D] = min(abs(DATA.Ts0D-tstart));
...@@ -11,11 +11,11 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS,CODE) ...@@ -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_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE;
Qx_infty_std = std (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)]); % disp(['Q_x=',sprintf('%2.2e',Qx_infty_avg),'+-',sprintf('%2.2e',Qx_infty_std)]);
f_avg_z = squeeze(mean(DATA.PHI(:,:,:,:),3)); % f_avg_z = squeeze(mean(DATA.PHI(:,:,:,:),3));
[~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3))); % [~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3)));
ikzf = min([ikzf,DATA.Nky]); % ikzf = min([ikzf,DATA.Nky]);
Ns3D = numel(DATA.Ts3D); % Ns3D = numel(DATA.Ts3D);
[KX, KY] = meshgrid(DATA.kx, DATA.ky); % [KX, KY] = meshgrid(DATA.kx, DATA.ky);
%% error estimation %% error estimation
DT_ = (tend-tstart)/OPTIONS.NCUT; DT_ = (tend-tstart)/OPTIONS.NCUT;
Qx_ee = zeros(1,OPTIONS.NCUT); Qx_ee = zeros(1,OPTIONS.NCUT);
...@@ -27,36 +27,36 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS,CODE) ...@@ -27,36 +27,36 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS,CODE)
Qx_avg = mean(Qx_ee); Qx_avg = mean(Qx_ee);
Qx_err = std(Qx_ee); Qx_err = std(Qx_ee);
disp(['Q_avg=',sprintf('%2.2e',Qx_avg),'+-',sprintf('%2.2e',Qx_err)]); disp(['Q_avg=',sprintf('%2.2e',Qx_avg),'+-',sprintf('%2.2e',Qx_err)]);
%% computations % %% computations
%
% Compute zonal and non zonal energies % % Compute zonal and non zonal energies
E_Zmode_SK = zeros(1,Ns3D); % E_Zmode_SK = zeros(1,Ns3D);
E_NZmode_SK = zeros(1,Ns3D); % E_NZmode_SK = zeros(1,Ns3D);
for it = 1:numel(DATA.Ts3D) % 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_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))))); % E_NZmode_SK(it) = squeeze(sum(sum(((1+KX.^2+KY.^2).*abs(squeeze(f_avg_z(:,:,it))).^2.*(KY~=0)))));
end % end
% Compute thermodynamic entropy Eq.(5) Navarro et al. 2012 PoP % % Compute thermodynamic entropy Eq.(5) Navarro et al. 2012 PoP
% 1/2 sum_p sum_j Napj^2(k=0) (avg z) % % 1/2 sum_p sum_j Napj^2(k=0) (avg z)
Nipjz = sum(sum(sum(sum(conj(DATA.Nipj).*DATA.Nipj)))); % Nipjz = sum(sum(sum(sum(conj(DATA.Nipj).*DATA.Nipj))));
ff = trapz(DATA.z,Nipjz,5); % ff = trapz(DATA.z,Nipjz,5);
E_TE = 0.5*squeeze(ff); % E_TE = 0.5*squeeze(ff);
% Compute electrostatic energy % % Compute electrostatic energy
E_ES = zeros(size(DATA.Ts5D)); % E_ES = zeros(size(DATA.Ts5D));
bi = sqrt(KX.^2+KY.^2)*DATA.sigma(1)*sqrt(2*DATA.tau(1)); %argument of the kernel % bi = sqrt(KX.^2+KY.^2)*DATA.sigma(1)*sqrt(2*DATA.tau(1)); %argument of the kernel
for it5D = 1:numel(DATA.Ts5D) % for it5D = 1:numel(DATA.Ts5D)
[~,it3D] = min(abs(DATA.Ts3D-DATA.Ts5D(it5D))); % [~,it3D] = min(abs(DATA.Ts3D-DATA.Ts5D(it5D)));
for in = 1:DATA.Jmaxi % for in = 1:DATA.Jmaxi
Knphi = kernel(in-1,bi).*squeeze(trapz(DATA.z,DATA.PHI(:,:,:,it3D),3)); % 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)); % Ni0n_z= squeeze(trapz(DATA.z,DATA.Nipj(1,in,:,:,:,it5D),5));
E_ES(it5D) = 0.5*sum(sum(abs(conj(Knphi).*Ni0n_z))); % E_ES(it5D) = 0.5*sum(sum(abs(conj(Knphi).*Ni0n_z)));
end % end
end % end
%% Figure %% Figure
clr_ = lines(20); clr_ = lines(20);
mvm = @(x) movmean(x,OPTIONS.NMVA); mvm = @(x) movmean(x,OPTIONS.NMVA);
FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; %set(gcf, 'Position', [500, 1000, 1000, 600]) 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),'--',... plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'--',...
'color',clr_((DATA.Pmaxi-1)/2,:),... 'color',clr_((DATA.Pmaxi-1)/2,:),...
'DisplayName',['$\Gamma_x$ ',DATA.paramshort]); hold on; 'DisplayName',['$\Gamma_x$ ',DATA.paramshort]); hold on;
...@@ -80,16 +80,16 @@ mvm = @(x) movmean(x,OPTIONS.NMVA); ...@@ -80,16 +80,16 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
title({DATA.param_title,... title({DATA.param_title,...
['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$, Q^{\infty} = $',num2str(Qx_infty_avg)]}); ['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$, Q^{\infty} = $',num2str(Qx_infty_avg)]});
%% Free energy % %% Free energy
FIGURE.ax2 = subplot(3,1,2,'parent',FIGURE.fig); % FIGURE.ax2 = subplot(3,1,2,'parent',FIGURE.fig);
yyaxis left % yyaxis left
plot(DATA.Ts5D,E_TE,'DisplayName','$\epsilon_f$'); hold on; % plot(DATA.Ts5D,E_TE,'DisplayName','$\epsilon_f$'); hold on;
ylabel('Entropy');%('$\epsilon_f$') % ylabel('Entropy');%('$\epsilon_f$')
yyaxis right % yyaxis right
plot(DATA.Ts5D,E_ES,'DisplayName','$\epsilon_\phi$'); % plot(DATA.Ts5D,E_ES,'DisplayName','$\epsilon_\phi$');
ylabel('ES energy');%('$\epsilon_\phi$') % ylabel('ES energy');%('$\epsilon_\phi$')
xlim([DATA.Ts5D(1), DATA.Ts5D(end)]); % xlim([DATA.Ts5D(1), DATA.Ts5D(end)]);
xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]); % xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
%% radial shear radial profile %% radial shear radial profile
% computation % computation
Ns3D = numel(DATA.Ts3D); Ns3D = numel(DATA.Ts3D);
...@@ -108,7 +108,7 @@ mvm = @(x) movmean(x,OPTIONS.NMVA); ...@@ -108,7 +108,7 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
f2plot = toplot.FIELD; f2plot = toplot.FIELD;
dframe = ite3D - its3D; dframe = ite3D - its3D;
clim = max(max(max(abs(plt(f2plot(:,:,:)))))); 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)); [TY,TX] = meshgrid(DATA.x,DATA.Ts3D(toplot.FRAMES));
pclr = pcolor(TX,TY,squeeze(plt(f2plot))'); pclr = pcolor(TX,TY,squeeze(plt(f2plot))');
set(pclr, 'edgecolor','none'); set(pclr, 'edgecolor','none');
......
matlab/profiler.m
+ 1
− 1
View file @ 11ff59c4
...@@ -95,7 +95,7 @@ end ...@@ -95,7 +95,7 @@ end
p1 = area(xx_,yy_,'LineStyle','none'); p1 = area(xx_,yy_,'LineStyle','none');
for i = 1:N_T; p1(i).FaceColor = colors(i,:); 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} = 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; end;
legend(LEGEND,'Location','bestoutside'); legend(LEGEND,'Location','bestoutside');
% legend('Compute RHS','Adv. fields','ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Process', 'Missing') % legend('Compute RHS','Adv. fields','ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Process', 'Missing')
......
wk/analysis_gyacomo.m
+ 5
− 6
View file @ 11ff59c4
...@@ -80,7 +80,7 @@ options.NAME = '\phi'; ...@@ -80,7 +80,7 @@ options.NAME = '\phi';
% options.NAME = 'Q_x'; % options.NAME = 'Q_x';
% options.NAME = 'n_i'; % options.NAME = 'n_i';
% options.NAME = 'n_i-n_e'; % options.NAME = 'n_i-n_e';
options.PLAN = 'xy'; options.PLAN = 'kxky';
% options.NAME = 'f_i'; % options.NAME = 'f_i';
% options.PLAN = 'sx'; % options.PLAN = 'sx';
options.COMP = 'avg'; options.COMP = 'avg';
...@@ -105,16 +105,15 @@ options.NORMALIZE = 0; ...@@ -105,16 +105,15 @@ options.NORMALIZE = 0;
% options.NAME = '\omega_z'; % options.NAME = '\omega_z';
% options.NAME = 'T_i'; % options.NAME = 'T_i';
% options.NAME = 'n_i'; % options.NAME = 'n_i';
% options.NAME = '\phi^{NZ}'; options.NAME = '\phi^{NZ}';
options.NAME = 'N_i^{00}'; % options.NAME = 'N_i^{00}';
% options.NAME = 'N_i^{00}-N_e^{00}'; % options.NAME = 'N_i^{00}-N_e^{00}';
% options.NAME = 's_{Ex}'; % options.NAME = 's_{Ex}';
% options.NAME = 'Q_x'; % options.NAME = 'Q_x';
% options.NAME = 'k^2n_e'; % options.NAME = 'k^2n_e';
options.PLAN = 'kxky'; options.PLAN = 'xy';
options.COMP = 'avg'; options.COMP = 'avg';
% options.TIME = [49 50 51]; options.TIME = [80 200];
options.TIME = data.Ts3D(51:55);
options.RESOLUTION = 256; options.RESOLUTION = 256;
data.a = data.EPS * 2e3; data.a = data.EPS * 2e3;
......
wk/fast_analysis.m
+ 40
− 7
View file @ 11ff59c4
...@@ -4,17 +4,50 @@ PARTITION = '/misc/gyacomo23_outputs/'; ...@@ -4,17 +4,50 @@ PARTITION = '/misc/gyacomo23_outputs/';
% PARTITION = gyacomodir; % PARTITION = gyacomodir;
%% CBC %% 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,'/']; DATADIR = [PARTITION,resdir,'/'];
data = {}; 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');
%% options.TAVG_0 = 100;
figure; options.TAVG_1 = 1000;
plot(data.Ts0D,data.HFLUX_X); options.NCUT = 5; % Number of cuts for averaging and error estimation
\ No newline at end of file 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')
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment