Homec4science
Diffusion GYACOMO (OBSOLETE: LATEST VERSION ON gitlab.epfl.ch) (master) Browse

gyacomo/matlab/plot/plot_radial_transport_and_spacetime.m
5bb85d4b05f5master

Last Change
Raw File

plot_radial_transport_and_spacetime.m
View Options

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));
[~,ite0D] = min(abs(DATA.Ts0D-tend));
[~,its3D] = min(abs(DATA.Ts3D-tstart));
[~,ite3D] = min(abs(DATA.Ts3D-tend));
SCALE = DATA.scale;%(1/DATA.Nx/DATA.Ny)^2;
Gx_infty_avg = mean(DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
Gx_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
Qx_infty_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE;
Qx_infty_std = std (DATA.HFLUX_X(its0D:ite0D))*SCALE;
disp(['G_x=',sprintf('%2.2e',Gx_infty_avg),'+-',sprintf('%2.2e',Gx_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));
[~,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);
for i = 1:OPTIONS.NCUT
[~,its_] = min(abs(DATA.Ts0D - (tstart+(i-1)*DT_)));
[~,ite_] = min(abs(DATA.Ts0D - (tstart+ i *DT_)));
Qx_ee(i) = mean(DATA.HFLUX_X(its_:ite_))*SCALE;
end
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 Gamma from ifft matlab
Gx = zeros(DATA.Ny,DATA.Nx,numel(DATA.Ts3D));
% for it = 1:numel(DATA.Ts3D)
% for iz = 1:DATA.Nz
% Gx(:,:,it) = Gx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
% .*ifourier_GENE(DATA.DENS_I(:,:,iz,it));
% end
% Gx(:,:,it) = Gx(:,:,it)/DATA.Nz;
% end
Gx_t_mtlb = squeeze(mean(mean(Gx,1),2));
% Compute Heat flux from ifft matlab
Qx = zeros(DATA.Ny,DATA.Nx,numel(DATA.Ts3D));
% for it = 1:numel(DATA.Ts3D)
% for iz = 1:DATA.Nz
% Qx(:,:,it) = Qx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
% .*ifourier_GENE(DATA.TEMP_I(:,:,iz,it));
% end
% Qx(:,:,it) = Qx(:,:,it)/DATA.Nz;
% end
Qx_t_mtlb = squeeze(mean(mean(Qx,1),2));
% zonal vs nonzonal energies for phi(t)
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
%% Figure
mvm = @(x) movmean(x,OPTIONS.NMVA);
FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [100, 100, 1000, 600])
subplot(311)
% yyaxis left
plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
% plot(mvm(DATA.Ts3D),mvm(Gx_t_mtlb),'DisplayName','matlab comp.'); hold on;
% plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Gx_infty_avg, '-k',...
% 'DisplayName',['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$\pm$',num2str(Gx_infty_std)]);
% ylabel('$\Gamma_x$')
% ylim([0,5*abs(Gx_infty_avg)]);
% xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
% yyaxis right
plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
% plot(mvm(DATA.Ts3D),mvm(Qx_t_mtlb),'DisplayName','matlab comp.'); hold on;
ylabel('Transport')
if(~isnan(Qx_infty_avg))
plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_infty_avg, '-k',...
'DisplayName',['$Q_{avg}=',sprintf('%2.2f',Qx_avg),'\pm',sprintf('%2.2f',Qx_err),'$']); legend('show');
ylim([0,5*abs(Qx_infty_avg)]);
else
plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Gx_infty_avg, '-k',...
'DisplayName',['$\Gamma_{avg}=',sprintf('%2.2f',Gx_infty_avg),'\pm',sprintf('%2.2f',Gx_infty_std),'$']); legend('show');
try ylim([0,5*abs(Gx_infty_avg)]);
catch
end
end
xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
grid on; set(gca,'xticklabel',[]);
title({DATA.param_title,...
['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$, Q^{\infty} = $',num2str(Qx_infty_avg)]});
%% radial shear radial profile
% computation
Ns3D = numel(DATA.Ts3D);
[KX, KY] = meshgrid(DATA.kx, DATA.ky);
plt = @(x) mean(x(:,:,:),1);
kycut = max(DATA.ky);
kxcut = max(DATA.kx);
LP = (abs(KY)<kycut).*(abs(KX)<kxcut); %Low pass filter
OPTIONS.NAME = OPTIONS.ST_FIELD;
OPTIONS.PLAN = 'xy';
OPTIONS.COMP = 'avg';
OPTIONS.TIME = DATA.Ts3D;
OPTIONS.POLARPLOT = 0;
toplot = process_field(DATA,OPTIONS);
f2plot = toplot.FIELD;
dframe = ite3D - its3D;
clim = max(max(max(abs(plt(f2plot(:,:,:))))));
subplot(313)
[TY,TX] = meshgrid(DATA.x,DATA.Ts3D(toplot.FRAMES));
pclr = pcolor(TX,TY,squeeze(plt(f2plot))');
set(pclr, 'edgecolor','none');
legend(['$\langle ',OPTIONS.ST_FIELD,' \rangle_{y,z}$']) %colorbar;
caxis(clim*[-1 1]);
cmap = bluewhitered(256);
colormap(cmap)
xlabel('$t c_s/R$'), ylabel('$x/\rho_s$');
if OPTIONS.INTERP
shading interp
end
if strcmp(OPTIONS.ST_FIELD,'Gx')
subplot(311)
plot(DATA.Ts3D,squeeze(mean(plt(f2plot),1)));
end
%% Zonal vs NZonal energies
subplot(312)
it0 = 1; itend = Ns3D;
trange = toplot.FRAMES;
plt1 = @(x) x;%-x(1);
plt2 = @(x) x./max((x(:)));
toplot = sum(squeeze(plt(f2plot)).^2,1); %ST from before
% plty = @(x) x(500:end)./max(squeeze(x(500:end)));
yyaxis left
% plot(plt1(DATA.Ts3D(trange)),plt2(E_Zmode_SK(trange)),'DisplayName','$k_{zf}^2|\phi_{kzf}|^2$');
plot(plt1(DATA.Ts3D(trange)),plt2(toplot(:)),'DisplayName','Sum $A^2$');
ylim([-0.1, 1.5]); ylabel('$E_{Z}$')
yyaxis right
plot(plt1(DATA.Ts3D(trange)),plt2(E_NZmode_SK(trange)),'DisplayName','$(1+k^2)|\phi^k|^2$');
xlim([DATA.Ts3D(it0), DATA.Ts3D(itend)]);
ylim([-0.1, 1.5]); ylabel('$E_{NZ}$')
xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
end
c4science · Help