diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 17b9a01892ae518e96bd9e84c8addadd131b0d49..f6ccd44ba8f871c5005f278bad5746eaf993885f 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -45,10 +45,18 @@ DATA.outfilenames = {};
ii = 1;
while(CONTINUE)
filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM);
+ % Check presence and jobnummax
if (exist(filename, 'file') == 2 && JOBNUM <= JOBNUMMAX)
DATA.outfilenames{ii} = filename;
+ %test if it is corrupted or currently running
+ try
+ openable = ~isempty(h5read(filename,'/data/var3d/time'));
+ catch
+ openable = 0
+ end
+ if openable
%% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- disp(['Loading ',filename])
+ disp(sprintf('Loading ID %.2d (%s)',JOBNUM,filename));
% Loading from output file
CPUTIME = h5readatt(filename,'/data/input','cpu_time');
DT_SIM = h5readatt(filename,'/data/input','dt');
@@ -219,6 +227,7 @@ while(CONTINUE)
LASTJOB = JOBNUM;
Pe_old = Pe_new; Je_old = Je_new;
Pi_old = Pi_new; Ji_old = Ji_new;
+ end
elseif (JOBNUM > JOBNUMMAX)
CONTINUE = 0;
disp(['found ',num2str(JOBFOUND),' results']);
diff --git a/matlab/compute/compute_3D_zpinch_growth_rate.m b/matlab/compute/compute_3D_zpinch_growth_rate.m
index c89597b199ae13fc9730969a80577c60710c94ca..76340d0737f32c4025f8c8c59814c644f84a114c 100644
--- a/matlab/compute/compute_3D_zpinch_growth_rate.m
+++ b/matlab/compute/compute_3D_zpinch_growth_rate.m
@@ -4,7 +4,7 @@ FIGURE.FIGNAME = ['growth_rate_kx=0_ky=0_planes',DATA.PARAMS];
t = DATA.Ts3D;
[~,its] = min(abs(t-TRANGE(1)));
[~,ite] = min(abs(t-TRANGE(end)));
-nkx = DATA.Nkx; nky = DATA.Nky; nkz = DATA.Nz/2;
+nkx = DATA.Nkx; nky = DATA.Nky; nt = numel(t);
% Remove AA part
if DATA.Nx > 1
ikxnz = abs(DATA.PHI(1,:,1,1)) > 0;
@@ -15,19 +15,27 @@ ikynz = abs(DATA.PHI(:,1,1,1)) > 0;
ikynz(1) = 1; ikxnz(1) = 1; %put k=0 in the analysis
+% phi = fft(DATA.PHI(:,:,:,:),[],3);
Y = fft(DATA.PHI(ikynz,ikxnz,:,:),[],3);
-phi = Y(:,:,2:2:end,:);
+phi = Y(:,:,2:2:end,:); sz_=size(phi);
+nkz = sz_(3);
+% tmp_ = ifourier_GENE(DATA.PHI(:,:,:,1)); sz_ = size(tmp_);
+% phi = zeros(sz_(1),sz_(2),sz_(3),nt);
+% for it = 1:nt
+% tmp_ = ifourier_GENE(DATA.PHI(:,:,:,it));
+% phi(:,:,:,it) = fftn(tmp_);
+% end
omega = zeros(nky,nkx,nkz);
for ikz = 1:nkz
for iky = 1:nky
- for ix = 1:nkx
+ for ikx = 1:nkx
% omega(iy,ix,iz) = LinearFit_s(t(its:ite),squeeze(abs(phi(iy,ix,iz,its:ite))));
- to_measure = squeeze(log(phi(iky,ix,ikz,its:ite)));
+ to_measure = squeeze(log(abs(phi(iky,ikx,ikz,its:ite))));
tmp = polyfit(t(its:ite),to_measure(:),1);
if ~(isnan(tmp(1)) || isinf(tmp(1)))
- omega(iky,ix,ikz) = tmp(1);
+ omega(iky,ikx,ikz) = tmp(1);
end
end
end
diff --git a/matlab/compute/compute_fa_1D.m b/matlab/compute/compute_fa_1D.m
index 1d24967649aaf29c14351c402520a18b2996613e..f8d569461615739f1bb93a576a53e4ebeef19288 100644
--- a/matlab/compute/compute_fa_1D.m
+++ b/matlab/compute/compute_fa_1D.m
@@ -17,7 +17,7 @@ xmin = min(abs(x));
[~, iky0] = min(abs(data.ky));
kx_ = data.kx; ky_ = data.ky;
-switch options.SPECIE
+switch options.SPECIES
case 'e'
Napj_ = data.Nepj;
parray = double(data.Pe);
diff --git a/matlab/compute/compute_fa_2D.m b/matlab/compute/compute_fa_2D.m
index a9c099318a1e3ccf7e3a17b46009e592e5beb5ba..8d9758e0b27a106faad7bce9a9bce4058ee67f57 100644
--- a/matlab/compute/compute_fa_2D.m
+++ b/matlab/compute/compute_fa_2D.m
@@ -15,7 +15,7 @@ FaM = @(s,x) exp(-s.^2-x);
[~, iky0] = min(abs(data.ky));
kx_ = data.kx; ky_ = data.ky;
-switch options.SPECIE
+switch options.SPECIES
case 'e'
Napj_ = data.Nepj;
parray = double(data.Pe);
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index dda3387ecda1f1598d7d0b95dc52f4fc88ccc3aa..076058a77e13513ab248ecdd88dbd45adff55642 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -2,12 +2,16 @@ function [FIGURE] = mode_growth_meter(DATA,OPTIONS)
NORMALIZED = OPTIONS.NORMALIZED;
Nma = OPTIONS.NMA; %Number moving average
+d = OPTIONS.fftz.flag; % To add spectral evolution of z (useful for 3d zpinch)
+
switch OPTIONS.iz
case 'avg'
field = squeeze(mean(DATA.PHI,3));
+ zstrcomp = 'z-avg';
otherwise
field = squeeze(DATA.PHI(:,:,OPTIONS.iz,:));
+ zstrcomp = ['z=',DATA.z(OPTIONS.iz)];
end
FRAMES = zeros(size(OPTIONS.TIME));
@@ -18,16 +22,18 @@ FRAMES = unique(FRAMES);
t = DATA.Ts3D(FRAMES);
% time window where we measure the growth
-it1 = floor(numel(t)/2);
-it2 = numel(t);
+TW = [OPTIONS.KY_TW; OPTIONS.KX_TW];
[~,ikzf] = max(mean(squeeze(abs(field(1,:,FRAMES))),2));
FIGURE.fig = figure; set(gcf, 'Position', [100 100 1200 700])
FIGURE.FIGNAME = 'mode_growth_meter';
for i = 1:2
- MODES_SELECTOR = i; %(2:Zonal, 1: NZonal)
+ MODES_SELECTOR = i; %(1:kx=0; 2:ky=0)
+ [~,it1] = min(abs(t-TW(i,1)));
+ [~,it2] = min(abs(t-TW(i,2)));
+
if MODES_SELECTOR == 2
switch OPTIONS.ik
case 'sum'
@@ -68,13 +74,6 @@ for i = 1:2
end
MODES = 1:numel(k);
- % MODES = zeros(size(OPTIONS.K2PLOT));
- % for i = 1:numel(OPTIONS.K2PLOT)
- % [~,MODES(i)] =min(abs(OPTIONS.K2PLOT(i)-k));
- % end
-
-
- % plt = @(x,ik) abs(squeeze(x(1,ik,iz,FRAMES)))./max(abs(squeeze(x(1,ik,iz,FRAMES))));
gamma = MODES;
amp = MODES;
@@ -89,15 +88,17 @@ for i = 1:2
mod2plot = [2:min(Nmax,OPTIONS.NMODES+1)];
clr_ = jet(numel(mod2plot));
%plot
- subplot(2,3,1+3*(i-1))
+% subplot(2+d,3,1+3*(i-1))
+ FIGURE.axes(1+3*(i-1)) = subplot(2+d,3,1+3*(i-1),'parent',FIGURE.fig);
[YY,XX] = meshgrid(t,fftshift(k));
pclr = pcolor(XX,YY,abs(plt(fftshift(field,MODES_SELECTOR),1:numel(k))));set(pclr, 'edgecolor','none'); hold on;
set(gca,'YDir','normal')
% xlim([t(1) t(end)]); %ylim([1e-5 1])
xlabel(['$',kstr,'\rho_s$']); ylabel('$t c_s /\rho_s$');
- title(MODESTR)
+ title([MODESTR,', ',zstrcomp])
- subplot(2,3,2+3*(i-1))
+% subplot(2+d,3,2+3*(i-1))
+ FIGURE.axes(2+3*(i-1)) = subplot(2+d,3,2+3*(i-1),'parent',FIGURE.fig);
for i_ = 1:numel(mod2plot)
semilogy(t,plt(field,MODES(mod2plot(i_))),'color',clr_(i_,:)); hold on;
% semilogy(t,exp(gamma(i_).*t+amp(i_)),'--k')
@@ -112,7 +113,8 @@ for i = 1:2
xlabel('$t c_s /\rho_s$'); ylabel(['$|\phi_{',kstr,'}|$']);
title('Measure time window')
- subplot(2,3,3+3*(i-1))
+% subplot(2+d,3,3+3*(i-1))
+ FIGURE.axes(3+3*(i-1)) = subplot(2+d,3,3+3*(i-1),'parent',FIGURE.fig);
plot(k(MODES),gamma,...
'DisplayName',['(',num2str(DATA.Pmaxi-1),',',num2str(DATA.Jmaxi-1),')']); hold on;
for i_ = 1:numel(mod2plot)
@@ -124,4 +126,60 @@ for i = 1:2
xlabel(['$',kstr,'\rho_s$']); ylabel('$\gamma$');
title('Growth rates')
end
+
+if d
+ [~,ikx] = min(abs(DATA.kx-OPTIONS.fftz.kx));
+ [~,iky] = min(abs(DATA.ky-OPTIONS.fftz.ky));
+ sz_=size(DATA.PHI);nkz = sz_(3)/2;
+ k = [(0:nkz/2), (-nkz/2+1):-1]/DATA.Npol;
+ % Spectral treatment of z-dimension
+ Y = fft(DATA.PHI(iky,ikx,:,:),[],3);
+ phi = squeeze(Y(1,1,2:2:end,:));
+
+ gamma = zeros(nkz);
+ for ikz = 1:nkz
+ to_measure = squeeze(log(abs(phi(ikz,it1:it2))));
+ tmp = polyfit(t(it1:it2),to_measure(:),1);
+ if ~(isnan(tmp(1)) || isinf(tmp(1)))
+ gamma(ikz) = tmp(1);
+ end
+ end
+ %plot
+ MODES = 1:numel(k);
+ mod2plot = [2:2:min(nkz,OPTIONS.NMODES+1)];
+ clr_ = jet(numel(mod2plot));
+ subplot(3,3,7)
+ [YY,XX] = meshgrid(t,fftshift(k));
+ pclr = pcolor(XX,YY,abs(phi));set(pclr, 'edgecolor','none'); hold on;
+ set(gca,'YDir','normal')
+ % xlim([t(1) t(end)]); %ylim([1e-5 1])
+ xlabel('$k_\parallel$'); ylabel('$t c_s /\rho_s$');
+ title(['$k_x=$',num2str(DATA.kx(ikx)),', $k_y=$',num2str(DATA.ky(iky))]);
+
+ subplot(3,3,8)
+ for i_ = 1:numel(mod2plot)
+ im_ = MODES(mod2plot(i_));
+ semilogy(t,abs(phi(im_,:)),'color',clr_(i_,:)); hold on;
+ end
+ %plot the time window where the gr are measured
+ plot(t(it1)*[1 1],[1e-10 1],'--k')
+ plot(t(it2)*[1 1],[1e-10 1],'--k')
+ xlim([t(1) t(end)]); %ylim([1e-5 1])
+ xlabel('$t c_s /\rho_s$'); ylabel(['$|\phi_{',kstr,'}|$']);
+ title('Measure time window')
+
+ subplot(3,3,9)
+ plot(k(MODES),gamma,...
+ 'DisplayName',['(',num2str(DATA.Pmaxi-1),',',num2str(DATA.Jmaxi-1),')']); hold on;
+ for i_ = 1:numel(mod2plot)
+ plot(k(MODES(mod2plot(i_))),gamma(mod2plot(i_)),'x','color',clr_(i_,:));
+ end
+ if MODES_SELECTOR == 2
+ plot(k(ikzf),gamma(ikzf),'ok');
+ end
+ xlabel(['$k_\parallel$']); ylabel('$\gamma$');
+ title('Growth rates')
+
+end
+
end
\ No newline at end of file
diff --git a/matlab/compute/process.m b/matlab/compute/process.m
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index 705b7f3a80de9d067a73494cdc31b30dd4c7405c..4be8de84db0369770847a6d81944549fb9433c44 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 = [00 40000];
+tw = [400 40000];
fig = gcf;
axObjs = fig.Children;
diff --git a/matlab/load/load_gene_data.m b/matlab/load/load_gene_data.m
index 641e6dd9aff55ea85d377d20fc3c1e984c3a9e22..63554d8a6647f2bd8c4669752f8a1875c6ecc3cc 100644
--- a/matlab/load/load_gene_data.m
+++ b/matlab/load/load_gene_data.m
@@ -2,6 +2,7 @@ function [ DATA ] = load_gene_data( folder )
%to load gene data as for HeLaZ results
namelist = read_namelist([folder,'parameters.dat']);
DATA.namelist = namelist;
+DATA.folder = folder;
%% Grid
coofile = 'coord.dat.h5';
DATA.vp = h5read([folder,coofile],'/coord/vp');
@@ -52,16 +53,24 @@ momfile = 'mom_ions.dat.h5';
for jt = 1:numel(DATA.Ts3D)
t = DATA.Ts3D(jt);
[~, it] = min(abs(DATA.Ts3D-t));
-%
-% tmp = h5read([folder,momfile],['/mom_ions/dens/',sprintf('%10.10d',it-1)]);
-% DATA.DENS_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
-% %
-% tmp = h5read([folder,momfile],['/mom_ions/T_par/',sprintf('%10.10d',it-1)]);
-% DATA.TPAR_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
-% %
-% tmp = h5read([folder,momfile],['/mom_ions/T_perp/',sprintf('%10.10d',it-1)]);
-% DATA.TPER_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
-%
+try
+ tmp = h5read([folder,momfile],['/mom_ions/dens/',sprintf('%10.10d',it-1)]);
+ DATA.DENS_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
+catch
+ DATA.DENS_I(:,:,:,it) = 0;
+end
+try
+ tmp = h5read([folder,momfile],['/mom_ions/T_par/',sprintf('%10.10d',it-1)]);
+ DATA.TPAR_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
+catch
+ DATA.TPAR_I(:,:,:,it) = 0;
+end
+try
+ tmp = h5read([folder,momfile],['/mom_ions/T_perp/',sprintf('%10.10d',it-1)]);
+ DATA.TPER_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
+catch
+ DATA.TPER_I(:,:,:,it) = 0;
+end
tmp = h5read([folder,phifile],['/field/phi/',sprintf('%10.10d',it-1)]);
DATA.PHI(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
diff --git a/matlab/plot/imagesc_custom.m b/matlab/plot/imagesc_custom.m
new file mode 100644
index 0000000000000000000000000000000000000000..3b57744e91760e2c17592b2ea526d358ceb38b9a
--- /dev/null
+++ b/matlab/plot/imagesc_custom.m
@@ -0,0 +1,12 @@
+function [ fig ] = imagesc_custom( XX,YY,FF)
+% CUSTOM IMAGESC this custom function is meant to be used exaclty as pcolor
+% but plotting array values in pixels (like imagesc) and not in vertices (like pcolor).
+
+fig = imagesc( XX(1,:), YY(:,1), FF);
+set(gca,'YDir','normal')
+xticks(XX(1,:));
+% xticklabels(num2str(XX(1,:)));
+yticks(YY(:,1));
+% yticklabels(num2str(YY(:,1)));
+end
+
diff --git a/matlab/plot/plot_fa.m b/matlab/plot/plot_fa.m
index 2395489c00ae1b8ddc38d69c99c60537ff193ae9..7faeabeba38cd7a9802a46badcaa2d0dad63bbde 100644
--- a/matlab/plot/plot_fa.m
+++ b/matlab/plot/plot_fa.m
@@ -10,14 +10,14 @@ end
if OPTIONS.ONED
[s,x,fsa,fxa] = compute_fa_1D(DATA, OPTIONS);
[~,it] = min(abs(OPTIONS.T-DATA.Ts5D));
- subplot(1,2,1)
+ FIGURE.ax1 = subplot(1,2,1,'parent',FIGURE.fig);
plot(s,fsa); hold on
- legend(OPTIONS.SPECIE)
+ legend(OPTIONS.SPECIES)
xlabel('$v_\parallel, (\mu=0)$'); ylabel(['$\langle |f_a|^2\rangle_{xy}^{1/2}$, ',zcomp]);
title(DATA.param_title);
- subplot(1,2,2)
+ FIGURE.ax2 = subplot(1,2,2,'parent',FIGURE.fig);
plot(x,fxa); hold on;
- legend(OPTIONS.SPECIE)
+ legend(OPTIONS.SPECIES)
xlabel('$\mu, (v_\parallel=0)$'); ylabel(['$\langle |f_a|^2\rangle_{xy}^{1/2}$, ',zcomp]);
if numel(it) == 1
title(['t=',num2str(DATA.Ts5D(it))]);
@@ -26,6 +26,7 @@ if OPTIONS.ONED
end
else
+ FIGURE.ax1 = subplot(1,1,1,'parent',FIGURE.fig);
[SS,XX,FFa] = compute_fa_2D(DATA, OPTIONS); sz = size(SS); FFa = FFa';
[~,it] = min(abs(OPTIONS.T-DATA.Ts5D));
switch OPTIONS.PLT_FCT
@@ -47,13 +48,13 @@ else
xlim([min(OPTIONS.SPAR) max(OPTIONS.SPAR)]);
ylim(sqrt(max(OPTIONS.XPERP))*[-1 1]);
end
- legend(['$\langle |f_',OPTIONS.SPECIE,'|^2\rangle_{xy}^{1/2}$',zcomp])
+ legend(['$\langle |f_',OPTIONS.SPECIES,'|^2\rangle_{xy}^{1/2}$',zcomp])
if numel(it) == 1
- title(['HeLaZ''$\langle |f_',OPTIONS.SPECIE...
+ title(['HeLaZ''$\langle |f_',OPTIONS.SPECIES...
,'|^2\rangle_{xy}^{1/2}$',zcomp...
,', $t=$[',sprintf('%1.1f',DATA.Ts5D(it)),']']);
else
- title(['HeLaZ''$\langle |f_',OPTIONS.SPECIE,...
+ title(['HeLaZ''$\langle |f_',OPTIONS.SPECIES,...
'|^2\rangle_{xy}^{1/2}$',zcomp,...
', average $t\in$[',sprintf('%1.1f',DATA.Ts5D(min(it))),',',sprintf('%1.1f',DATA.Ts5D(max(it))),']']);
end
diff --git a/matlab/plot/plot_fa_gene.m b/matlab/plot/plot_fa_gene.m
index 9743d80a3615d5b8e2a38239985c97851069b588..edcc073812d1f10b7d2edc2e1ae1fbca171d481d 100644
--- a/matlab/plot/plot_fa_gene.m
+++ b/matlab/plot/plot_fa_gene.m
@@ -1,7 +1,7 @@
-function plot_fa_gene(OPTIONS)
+function [FIGURE] = plot_fa_gene(OPTIONS)
folder = OPTIONS.folder;
-TIMES = OPTIONS.times;
-specie = OPTIONS.specie;
+TIMES = OPTIONS.T;
+SPECIES = OPTIONS.SPECIES;
PLT_FCT= OPTIONS.PLT_FCT;
file = 'coord.dat.h5';
@@ -25,7 +25,7 @@ end
file = 'vsp.dat.h5';
time = h5read([folder,file],'/vsp/time');
-fig = figure;
+FIGURE.fig = figure;
G_t = [];
Gdata = 0;
@@ -38,30 +38,32 @@ end
Gdata = Gdata ./ numel(TIMES);
if OPTIONS.ONED
- switch specie
+ switch SPECIES
case 'e'
FFa = zcomp(Gdata(:,:,:,2));
FFa = abs(FFa)./max(max(abs(FFa)));
- subplot(1,2,1)
+ FIGURE.ax1 = subplot(1,2,1,'parent',FIGURE.fig);
plot(vp,FFa(:,im0)); hold on;
- subplot(1,2,2)
+ FIGURE.ax2 = subplot(1,2,2,'parent',FIGURE.fig);
plot(mu,FFa(iv0,:)); hold on;
case 'i'
FFa = zcomp(Gdata(:,:,:,1));
FFa = abs(FFa)./max(max(abs(FFa)));
end
- subplot(1,2,1)
+ FIGURE.ax1 = subplot(1,2,1,'parent',FIGURE.fig);
plot(vp,FFa(:,im0)); hold on;
- legend(specie)
+ legend(SPECIES)
xlabel('$v_\parallel, (\mu=0)$'); ylabel('$\langle |f_a|^2\rangle_{xy}^{1/2}$');
- subplot(1,2,2)
+ FIGURE.ax2 = subplot(1,2,2,'parent',FIGURE.fig);
plot(mu,FFa(iv0,:)); hold on;
- legend(specie)
- xlabel('$\mu, (v_\parallel=0)$'); ylabel('$\langle |f_a|^2\rangle_{xy}^{1/2}$');
+ legend(SPECIES)
+ xlabel('$\mu, (v_\parallel=0)$'); %ylabel('$\langle |f_a|^2\rangle_{xy}^{1/2}$');
+ xticks(FIGURE.ax2,[]);
else
-switch specie
+ FIGURE.ax1 = subplot(1,1,1,'parent',FIGURE.fig);
+switch SPECIES
case 'e'
name = '$f_e(v_\parallel,\mu_p)$';
FFa = zcomp(squeeze(Gdata(:,:,:,2)));
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index ca11c1fd778ebf6406932882f3054f9bd35033ac..e37e1073dd2761be260846e32d2c2c4bcfba17b4 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)
+function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS,CODE)
%Compute steady radial transport
tend = OPTIONS.TAVG_1; tstart = OPTIONS.TAVG_0;
[~,its0D] = min(abs(DATA.Ts0D-tstart));
@@ -31,51 +31,41 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
% 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)
-
+ % 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)
+ switch CODE
+ case 'GYACOMO'
+ 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_i*sqrt(2*DATA.tau_i); %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
+ otherwise
+ E_TE = 0; E_ES =0; DATA.Ts5D =[0 1];
+ end
%% Figure
mvm = @(x) movmean(x,OPTIONS.NMVA);
- FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [500, 1000, 1000, 600])
- subplot(311)
-% yyaxis left
+ 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);
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',...
@@ -92,15 +82,26 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
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
+
+ %% 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);
[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';
@@ -110,7 +111,7 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
f2plot = toplot.FIELD;
dframe = ite3D - its3D;
clim = max(max(max(abs(plt(f2plot(:,:,:))))));
- subplot(313)
+ FIGURE.ax3 = subplot(3,1,3,'parent',FIGURE.fig);
[TY,TX] = meshgrid(DATA.x,DATA.Ts3D(toplot.FRAMES));
pclr = pcolor(TX,TY,squeeze(plt(f2plot))');
set(pclr, 'edgecolor','none');
@@ -126,21 +127,21 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
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]);
+ %% 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
\ No newline at end of file
diff --git a/matlab/plot/spectrum_1D.m b/matlab/plot/spectrum_1D.m
index c0669f239ce4bb72a9a32e4cad89cf1b7375b0e3..0bf0d8e1e127e58aecb70f7be73c6ae79482c708 100644
--- a/matlab/plot/spectrum_1D.m
+++ b/matlab/plot/spectrum_1D.m
@@ -66,8 +66,7 @@ end
if ~PLOT2D
set(gcf, 'Position', [20 50 1200 500])
- subplot(1,2,1)
-
+ FIGURE.ax1 = subplot(1,2,1,'parent',FIGURE.fig);
k = data.ky;
xname = '$k_y$';
@@ -111,7 +110,7 @@ if ~PLOT2D
% legend('show','Location','eastoutside')
xlabel(xname); ylabel(ynamecx)
- subplot(1,2,2)
+ FIGURE.ax2 = subplot(1,2,2,'parent',FIGURE.fig);
k = data.kx;
xname = '$k_x$';
diff --git a/matlab/process_field_v2.m b/matlab/process_field_v2.m
index 48122855a277532ef8f9aac340e42cd1f56dfa65..f744d7e49b760d42bbcaedee4ec31dd00ac3e1d5 100644
--- a/matlab/process_field_v2.m
+++ b/matlab/process_field_v2.m
@@ -10,7 +10,7 @@ end
FRAMES = unique(FRAMES);
%% Setup the plot geometry
[KX, KY] = meshgrid(DATA.kx, DATA.ky);
-directions = {'x','y','z'};
+directions = {'y','x','z'};
Nx = DATA.Nx; Ny = DATA.Ny; Nz = DATA.Nz; Nt = numel(FRAMES);
POLARPLOT = OPTIONS.POLARPLOT;
LTXNAME = OPTIONS.NAME;
@@ -39,6 +39,14 @@ switch OPTIONS.PLAN
XNAME = '$k_y$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.ky);
REALP = 0; COMPDIM = 2; POLARPLOT = 0; SCALE = 0;
+ case 'sx'
+ XNAME = '$v_\parallel$'; YNAME = '$\mu$';
+ [Y,X] = meshgrid(OPTIONS.XPERP,OPTIONS.SPAR);
+ REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 0;
+ for i = 1:numel(OPTIONS.TIME)
+ [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts5D));
+ end
+ FRAMES = unique(FRAMES); Nt = numel(FRAMES);
end
Xmax_ = max(max(abs(X))); Ymax_ = max(max(abs(Y)));
Lmin_ = min([Xmax_,Ymax_]);
@@ -66,27 +74,28 @@ else
end
%% Process the field to plot
% short term writing --
-b_e = DATA.sigma_e*sqrt(2*DATA.tau_e)*sqrt(KX.^2+KY.^2);
-adiab_e = kernel(0,b_e);
-pol_e = kernel(0,b_e).^2;
-for n = 1:DATA.Jmaxe
- adiab_e = adiab_e + kernel(n,b_e).^2;
- pol_e = pol_e + kernel(n,b_e).^2;
-end
-adiab_e = DATA.q_e/DATA.tau_e .* adiab_e;
-pol_e = DATA.q_e^2/DATA.tau_e * (1 - pol_e);
-
-b_i = DATA.sigma_i*sqrt(2*DATA.tau_i)*sqrt(KX.^2+KY.^2);
-adiab_i = kernel(0,b_i);
-pol_i = kernel(0,b_i).^2;
-for n = 1:DATA.Jmaxi
- adiab_i = adiab_i + kernel(n,b_i).^2;
- pol_i = pol_i + kernel(n,b_i).^2;
-end
-pol_i = DATA.q_i^2/DATA.tau_i * (1 - pol_i);
-adiab_i = DATA.q_i/DATA.tau_i .* adiab_i;
-poisson_op = (pol_i + pol_e);
-
+% b_e = DATA.sigma_e*sqrt(2*DATA.tau_e)*sqrt(KX.^2+KY.^2);
+% adiab_e = kernel(0,b_e);
+% pol_e = kernel(0,b_e).^2;
+% for n = 1:DATA.Jmaxe
+% adiab_e = adiab_e + kernel(n,b_e).^2;
+% pol_e = pol_e + kernel(n,b_e).^2;
+% end
+% adiab_e = DATA.q_e/DATA.tau_e .* adiab_e;
+% pol_e = DATA.q_e^2/DATA.tau_e * (1 - pol_e);
+%
+% b_i = DATA.sigma_i*sqrt(2*DATA.tau_i)*sqrt(KX.^2+KY.^2);
+% adiab_i = kernel(0,b_i);
+% pol_i = kernel(0,b_i).^2;
+% for n = 1:DATA.Jmaxi
+% adiab_i = adiab_i + kernel(n,b_i).^2;
+% pol_i = pol_i + kernel(n,b_i).^2;
+% end
+% pol_i = DATA.q_i^2/DATA.tau_i * (1 - pol_i);
+% adiab_i = DATA.q_i/DATA.tau_i .* adiab_i;
+% poisson_op = (pol_i + pol_e);
+adiab_e =0; adiab_i =0; poisson_op=1;
+SKIP_COMP = 0; % Turned on only for kin. distr. func. plot
% --
switch OPTIONS.COMP
case 'avg'
@@ -247,11 +256,34 @@ switch OPTIONS.NAME
end
FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.Nky,1:DATA.Nkx,:)));
end
+ case 'f_i'
+ SKIP_COMP = 1;
+ shift_x = @(x) x; shift_y = @(y) y;
+ NAME = 'fi'; OPTIONS.SPECIE = 'i';
+ for it = 1:numel(FRAMES)
+ OPTIONS.T = DATA.Ts5D(FRAMES(it));
+ OPTIONS.Z = OPTIONS.COMP;
+ [~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
+ end
+ case 'f_e'
+ SKIP_COMP = 1;
+ shift_x = @(x) x; shift_y = @(y) y;
+ NAME = 'fe'; OPTIONS.SPECIE = 'e';
+ [~,it0_] =min(abs(OPTIONS.TIME(1)-DATA.Ts5D));
+ [~,it1_]=min(abs(OPTIONS.TIME(end)-DATA.Ts5D));
+ dit_ = 1;%ceil((it1_-it0_+1)/10);
+ FRAMES = it0_:dit_:it1_;
+ iz = 1;
+ for it = 1:numel(FRAMES)
+ OPTIONS.T = DATA.Ts5D(FRAMES(it));
+ [~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
+ end
otherwise
disp('Fieldname not recognized');
return
end
% Process the field according to the 2D plane and the space (real/cpx)
+if ~SKIP_COMP
if COMPDIM == 3
for it = 1:numel(FRAMES)
tmp = squeeze(compr(OPE_.*FLD_(:,:,:,it)));
@@ -265,11 +297,12 @@ else
end
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(OPE_.*FLD_(:,:,:,it)));
+ tmp(:,:,iz) = squeeze(process(OPE_.*FLD_(:,:,iz,it)));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
+end
TOPLOT.FIELD = FIELD;
TOPLOT.FRAMES = FRAMES;
TOPLOT.X = shift_x(X);