diff --git a/matlab/LinearFit_s.m b/matlab/LinearFit_s.m
index dd5f11290c489beb4f1e3ee22541729b1002160d..6712007afeb58c49b1fadc0bd92599930c49ca77 100644
--- a/matlab/LinearFit_s.m
+++ b/matlab/LinearFit_s.m
@@ -22,7 +22,7 @@ function [gamma,fit] = LinearFit_s(time,Na00abs)
Na00absshifted = circshift(Na00abs,-1); % ... shift by -1 the time position
- gammaoft = log(Na00absshifted(1:end-1)./Na00abs(1:end-1))./squeeze(diff(time')); % ... evaluate growth rate
+ gammaoft = log(Na00absshifted(1:end-1)./Na00abs(1:end-1))./squeeze(diff(time)); % ... evaluate growth rate
% Get gamma
% gamma = mean(gammaoft(Tstart_ind:Tend_ind-1)); % ... take the mean of gamma over the time window
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 2200bbb05eec22902a58a33da4e809c1c157b027..ecf557383e32e7c152ab778a4812ae57c3cc229c 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -4,8 +4,7 @@ TJOB_SE = []; % Start and end times of jobs
NU_EVOL = []; % evolution of parameter nu between jobs
CO_EVOL = []; % evolution of CO
MU_EVOL = []; % evolution of parameter mu between jobs
-ETAN_EVOL= []; %
-ETAB_EVOL= []; %
+K_N_EVOL = []; %
L_EVOL = []; %
DT_EVOL = []; %
% FIELDS
@@ -25,7 +24,7 @@ Ts5D_ = [];
Sipj_ = []; Sepj_ = [];
Pe_old = 1e9; Pi_old = Pe_old; Je_old = Pe_old; Ji_old = Pe_old;
Pi_max=0; Pe_max=0; Ji_max=0; Je_max=0;
-while(CONTINUE)
+while(CONTINUE)
filename = sprintf([BASIC.MISCDIR,'outputs_%.2d.h5'],JOBNUM);
if (exist(filename, 'file') == 2 && JOBNUM <= JOBNUMMAX)
% Load results of simulation #JOBNUM
@@ -72,22 +71,22 @@ while(CONTINUE)
tmp = Sepj; sz = size(tmp);
Sepj = zeros(cat(1,[Pe_max,Je_max]',sz(3:end)')');
Sepj(1:Pe_new,1:Je_new,:,:,:) = tmp;
- end
+ end
end
-
+
if W_GAMMA || W_HF
Ts0D_ = cat(1,Ts0D_,Ts0D);
- end
-
+ end
+
if W_GAMMA
GGAMMA_ = cat(1,GGAMMA_,GGAMMA_RI);
PGAMMA_ = cat(1,PGAMMA_,PGAMMA_RI);
end
-
+
if W_HF
HFLUX_ = cat(1,HFLUX_,HFLUX_X);
end
-
+
if W_PHI || W_NA00
Ts3D_ = cat(1,Ts3D_,Ts3D);
end
@@ -120,15 +119,14 @@ while(CONTINUE)
% Evolution of simulation parameters
load_params
- TJOB_SE = [TJOB_SE Ts0D(1) Ts0D(end)];
+ TJOB_SE = [TJOB_SE Ts0D(1) Ts0D(end)];
NU_EVOL = [NU_EVOL NU NU];
CO_EVOL = [CO_EVOL CO CO];
MU_EVOL = [MU_EVOL MU MU];
- ETAN_EVOL = [ETAN_EVOL ETAN ETAN];
- ETAB_EVOL = [ETAB_EVOL ETAB ETAB];
+ K_N_EVOL = [K_N_EVOL K_N K_N];
L_EVOL = [L_EVOL L L];
DT_EVOL = [DT_EVOL DT_SIM DT_SIM];
-
+
JOBFOUND = JOBFOUND + 1;
LASTJOB = JOBNUM;
Pe_old = Pe_new; Je_old = Je_new;
@@ -148,4 +146,4 @@ clear Nipj_ Nepj_ Ni00_ Ne00_ PHI_ Ts2D_ Ts5D_ GGAMMA_ PGAMMA_ Ts0D_ HFLUX_
Sipj = Sipj_; Sepj = Sepj_;
clear Sipj_ Sepj_
JOBNUM = LASTJOB;
-filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM);
\ No newline at end of file
+filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM);
diff --git a/matlab/create_mov.m b/matlab/create_mov.m
index ca28c3c387d6aa30966b62722a847009a0d2d9eb..5f84d65d62ef1731ae21db7a34f862af1294933d 100644
--- a/matlab/create_mov.m
+++ b/matlab/create_mov.m
@@ -1,27 +1,28 @@
title1 = GIFNAME;
FIGDIR = BASIC.RESDIR;
-GIFNAME = [FIGDIR, GIFNAME,'.mp4'];
+GIFNAME = [FIGDIR, GIFNAME];
XNAME = latexize(XNAME);
YNAME = latexize(YNAME);
FIELDNAME = latexize(FIELDNAME);
vidfile = VideoWriter(GIFNAME,'Uncompressed AVI');
+vidfile.FrameRate = FPS;
open(vidfile);
% Set colormap boundaries
-hmax = max(max(max(FIELD)));
-hmin = min(min(min(FIELD)));
-
+hmax = max(max(max(FIELD(:,:,FRAMES))));
+hmin = min(min(min(FIELD(:,:,FRAMES))));
+scale = -1;
flag = 0;
if hmax == hmin
disp('Warning : h = hmin = hmax = const')
else
% Setup figure frame
-figure('Color','white','Position', [100, 100, 400, 400]);
+figure('Color','white','Position', [100, 100, 500, 500]);
pcolor(X,Y,FIELD(:,:,1)); % to set up
% colormap gray
colormap(bluewhitered)
axis tight manual % this ensures that getframe() returns a consistent size
- if 1
+ if INTERP
shading interp;
end
in = 1;
@@ -29,11 +30,23 @@ figure('Color','white','Position', [100, 100, 400, 400]);
for n = FRAMES % loop over selected frames
scale = max(max(abs(FIELD(:,:,n))));
pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot
- if 1
+ if INTERP
shading interp;
end
+ if NORMALIZED
+ pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot\
+ caxis([-1,1]);
+ else
+ pclr = pcolor(X,Y,FIELD(:,:,n)); % frame plot\
+ if CONST_CMAP
+ caxis([-1,1]*max(abs([hmin hmax]))); % const color map
+ else
+ caxis([-1,1]*scale); % adaptive color map
+ end
+ title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
+ ,', scale = ',sprintf('%.1e',scale)]);
+ end
set(pclr, 'edgecolor','none'); axis square;
-% caxis([-max(max(abs(FIELD(:,:,n)))),max(max(abs(FIELD(:,:,n))))]);
xlabel(XNAME); ylabel(YNAME); %colorbar;
title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
,', scaling = ',sprintf('%.1e',scale)]);
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
new file mode 100644
index 0000000000000000000000000000000000000000..65b13c91548e2b648e6f6c5281d70b7765f53ce8
--- /dev/null
+++ b/matlab/extract_fig_data.m
@@ -0,0 +1,10 @@
+fig = gcf;
+axObjs = fig.Children;
+dataObjs = axObjs.Children;
+
+X_ = dataObjs(1).XData;
+Y_ = dataObjs(1).YData;
+
+figure;
+plot(X_,Y_);
+plot(X_(9000:12000)-X_(8000),Y_(9000:12000));
\ No newline at end of file
diff --git a/matlab/load_marconi.m b/matlab/load_marconi.m
index 10406e8d8fb754dfcba7c038df5bc01f8fcd7cb2..faee436920e8fe661e0289cb7919c0f2403b51f8 100644
--- a/matlab/load_marconi.m
+++ b/matlab/load_marconi.m
@@ -1,14 +1,20 @@
function [ RESDIR ] = load_marconi( outfilename )
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
- hostfolder = outfilename(1:end-8);
- hostfile = [hostfolder,'/out*'];
- MISCDIR = ['/misc/HeLaZ_outputs/',outfilename(46:end-8),'/'];
- RESDIR = ['../',outfilename(46:end-8),'/'];
+ hostfolder = outfilename;
+ hostfile = [hostfolder,'out*'];
+ MISCDIR = ['/misc/HeLaZ_outputs/',outfilename(46:end)];
+ RESDIR = ['../',outfilename(46:end)];
miscfolder = [MISCDIR,'.'];
system(['mkdir -p ',miscfolder]);
disp(['mkdir -p ',miscfolder]);
resultfolder = [RESDIR,'.'];
+% CMD = ['rsync -r ahoffman@login.marconi.cineca.it:',hostfolder,'/* ',MISCDIR];
+% disp(CMD);
+% system(CMD);
+% CMD = ['rsync -r --exclude ''outputs_*.h5'' ahoffman@login.marconi.cineca.it:',hostfolder,'/* ',MISCDIR];
+% disp(CMD);
+% system(CMD);
% SCP the output file from marconi to misc folder of SPCPC
CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfile,' ',miscfolder];
disp(CMD);
diff --git a/matlab/load_params.m b/matlab/load_params.m
index 1e93d9e76e25a517dda6ea5777f661c7a48ef3cc..1b0abf43828ebee86a1e80bd1bcf018a78873712 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -1,7 +1,9 @@
CO = h5readatt(filename,'/data/input','CO');
-ETAB = h5readatt(filename,'/data/input','eta_B');
-ETAN = h5readatt(filename,'/data/input','eta_n');
-ETAT = h5readatt(filename,'/data/input','eta_T');
+% K_N = h5readatt(filename,'/data/input','eta_n');
+% K_T = h5readatt(filename,'/data/input','eta_T');
+K_N = h5readatt(filename,'/data/input','K_n');
+K_T = h5readatt(filename,'/data/input','K_T');
+K_E = h5readatt(filename,'/data/input','K_E');
PMAXI = h5readatt(filename,'/data/input','pmaxi');
JMAXI = h5readatt(filename,'/data/input','jmaxi');
PMAXE = h5readatt(filename,'/data/input','pmaxe');
@@ -49,9 +51,9 @@ else
degngrad = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),...
'_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI)];
end
-degngrad = [degngrad,'_eta_%1.1f_nu_%0.0e_',...
+degngrad = [degngrad,'_Kn_%1.1f_nu_%0.0e_',...
CONAME,'_CLOS_',num2str(CLOS),'_mu_%0.0e'];
-degngrad = sprintf(degngrad,[ETAB/ETAN,NU,MU]);
+degngrad = sprintf(degngrad,[K_N,NU,MU]);
if ~NON_LIN; degngrad = ['lin_',degngrad]; end
resolution = [num2str(Nx),'x',num2str(Ny/2),'_'];
gridname = ['L_',num2str(L),'_'];
diff --git a/matlab/load_results.m b/matlab/load_results.m
index 4ccc0671856229edefdc121e6cab796f27b83904..81bac625891aec776f3deb0eb56b96651937bc16 100644
--- a/matlab/load_results.m
+++ b/matlab/load_results.m
@@ -7,7 +7,7 @@ DT_SIM = h5readatt(filename,'/data/input','dt');
[Pe, Je, Pi, Ji, kx, ky, z] = load_grid_data(filename);
W_GAMMA = strcmp(h5readatt(filename,'/data/input','write_gamma'),'y');
-W_HF = strcmp(h5readatt(filename,'/data/input','write_hf' ),'y');
+W_HF = 0;%strcmp(h5readatt(filename,'/data/input','write_hf' ),'y');
W_PHI = strcmp(h5readatt(filename,'/data/input','write_phi' ),'y');
W_NA00 = strcmp(h5readatt(filename,'/data/input','write_Na00' ),'y');
W_NAPJ = strcmp(h5readatt(filename,'/data/input','write_Napj' ),'y');
diff --git a/matlab/plot_param_evol.m b/matlab/plot_param_evol.m
index 718bad2b3319f008ef964f27565a372150ad00ea..827b4c931985bf1db60ba83071f215b392152840 100644
--- a/matlab/plot_param_evol.m
+++ b/matlab/plot_param_evol.m
@@ -3,37 +3,37 @@ hold on; set(gcf, 'Position', [100, 100, 400, 1500])
subplot(4,1,1);
title('Parameter evolution'); hold on;
yyaxis left
- plot(TJOB_SE,NU_EVOL,'DisplayName','$\nu$');
+ plot(TJOB_SE,NU_EVOL,'DisplayName','$\nu$');
yyaxis right
- plot(TJOB_SE,CO_EVOL,'DisplayName','CO');
+ plot(TJOB_SE,CO_EVOL,'DisplayName','CO');
xlim([TJOB_SE(1) TJOB_SE(end)]);legend('show');grid on;
xticks([]);
- plot_tjob_lines(TJOB_SE,ylim)
+ plot_tjob_lines(TJOB_SE,ylim)
subplot(4,1,2);
plot(TJOB_SE,MU_EVOL,'DisplayName','$\mu$'); hold on;
xlim([TJOB_SE(1) TJOB_SE(end)]);legend('show');grid on;
xticks([]);
- plot_tjob_lines(TJOB_SE,ylim)
-
+ plot_tjob_lines(TJOB_SE,ylim)
+
subplot(4,1,3);
yyaxis left
- plot(TJOB_SE,ETAN_EVOL,'DisplayName','$\nabla n$'); hold on;
+ plot(TJOB_SE,K_N_EVOL,'DisplayName','$\nabla n$'); hold on;
yyaxis right
- plot(TJOB_SE,ETAB_EVOL,'DisplayName','$\nabla B$'); hold on;
+ plot(TJOB_SE,K_T_EVOL,'DisplayName','$\nabla T$'); hold on;
xlim([TJOB_SE(1) TJOB_SE(end)]);legend('show');grid on;
xticks([]);
- plot_tjob_lines(TJOB_SE,ylim)
+ plot_tjob_lines(TJOB_SE,ylim)
subplot(4,1,4);
plot(TJOB_SE,L_EVOL,'DisplayName','$L$'); hold on;
xlim([TJOB_SE(1) TJOB_SE(end)]);legend('show');grid on;
- plot_tjob_lines(TJOB_SE,ylim)
+ plot_tjob_lines(TJOB_SE,ylim)
xlabel('$t c_s/R$');
saveas(fig,[BASIC.RESDIR,'param_evol.png']);
-
+
function [] = plot_tjob_lines(TJOB,limits)
for i = 2:numel(TJOB)-1
- plot(TJOB(i)*[1 1],limits,'--k')
+ plot(TJOB(i)*[1 1],limits,'--k')
end
end
diff --git a/matlab/plots/plot_kperp_spectrum.m b/matlab/plots/plot_kperp_spectrum.m
index 43eadb463b9df1ae0ee033e127b6afa081c8b30f..10be9920f02c20757c08c025b376b959f6730c6d 100644
--- a/matlab/plots/plot_kperp_spectrum.m
+++ b/matlab/plots/plot_kperp_spectrum.m
@@ -29,13 +29,13 @@ plot(kp_ip,a1*kp_ip.^(-13/3),'-','DisplayName','$k^{-13/3}$');
plot(kp_ip,a2/100*kp_ip.^(-3)./(1+kp_ip.^2).^2,'-','DisplayName','$k^{-3}/(1+k^2)^2$');
end
if LOGSCALE
- set(gca, 'XScale', 'log');set(gca, 'YScale', 'log');
+ set(gca, 'XScale', 'log');set(gca, 'YScale', 'log');
xlim([0.1,10]);
end
grid on
xlabel('$k_\perp \rho_s$'); legend('show','Location','northeast')
-title({['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
+title({['$\nu_{',CONAME,'}=$', num2str(NU), ', $\kappa_N=$',num2str(K_N),...
', $L=',num2str(L),'$, $N=',num2str(Nx),'$'];[' $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
' $\mu_{hd}=$',num2str(MU),', $',num2str(round(tstart)),'<t<',num2str(round(tend)),'$']});
save_figure
-clear Z_rth phi_kp ni_kp Ti_kp
\ No newline at end of file
+clear Z_rth phi_kp ni_kp Ti_kp
diff --git a/matlab/plots/plot_radial_transport_and_shear.m b/matlab/plots/plot_radial_transport_and_shear.m
index 86bd431911cf1e85eaf6bf1c80930ea015efb5de..a8cdddeb655d2fb421eacef08ce2154308dacfec 100644
--- a/matlab/plots/plot_radial_transport_and_shear.m
+++ b/matlab/plots/plot_radial_transport_and_shear.m
@@ -5,6 +5,9 @@ tend = TAVG_1; tstart = TAVG_0;
SCALE = (2*pi/Nx/Ny)^2;
gamma_infty_avg = mean(PGAMMA_RI(its0D:ite0D))*SCALE;
gamma_infty_std = std (PGAMMA_RI(its0D:ite0D))*SCALE;
+if W_HF
+ q_infty_avg = mean(HFLUX_X(its0D:ite0D))*SCALE;
+end
% Compute steady shearing rate
tend = TAVG_1; tstart = TAVG_0;
[~,its2D] = min(abs(Ts3D-tstart));
@@ -14,32 +17,57 @@ shear_infty_std = std (mean(shear_maxx_avgy(:,its2D:ite2D),1));
% plots
fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
subplot(311)
-% yyaxis left
+ yyaxis left
plot(Ts0D,PGAMMA_RI*SCALE,'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
plot(Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*gamma_infty_avg, '-k',...
'DisplayName',['$\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]);
grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_x$')
ylim([0,5*abs(gamma_infty_avg)]); xlim([Ts0D(1),Ts0D(end)]);
- title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
+ title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\kappa_N=$',num2str(K_N),...
', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
' $\mu_{hd}=$',num2str(MU),', $\Gamma^{\infty} \approx $',num2str(gamma_infty_avg)]);
- %
+ if W_HF
+ yyaxis right
+ plot(Ts0D,HFLUX_X*SCALE,'DisplayName','$\langle T_i \partial_y\phi \rangle_y$'); hold on;
+ grid on; set(gca,'xticklabel',[]); ylabel('$Q_x$')
+ xlim([Ts0D(1),Ts0D(end)]); ylim([0,5*abs(q_infty_avg)]);
+ end
+ %
+% subplot(312)
+% clr = line_colors(1,:);
+% lstyle = line_styles(1);
+% % plt = @(x_) mean(x_,1);
+% plt = @(x_) x_(1,:);
+% plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{x,y}(\partial^2_x\phi)$'); hold on;
+% plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{x}\langle \partial^2_x\phi\rangle_y$'); hold on;
+% plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{y}\langle \partial^2_x\phi\rangle_x$'); hold on;
+% plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle \partial^2_x\phi\rangle_{x,y}$'); hold on;
+% plot(Ts3D(its2D:ite2D),ones(ite2D-its2D+1,1)*shear_infty_avg, '-k',...
+% 'DisplayName',['$s^{\infty} = $',num2str(shear_infty_avg),'$\pm$',num2str(shear_infty_std)]);
+% ylim([0,shear_infty_avg*5.0]); xlim([Ts0D(1),Ts0D(end)]);
+% grid on; ylabel('Shear amp.');set(gca,'xticklabel',[]);% legend('show');
+ %% zonal vs nonzonal energies for phi(t)
subplot(312)
- clr = line_colors(1,:);
- lstyle = line_styles(1);
-% plt = @(x_) mean(x_,1);
- plt = @(x_) x_(1,:);
- plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{x,y}(\partial^2_x\phi)$'); hold on;
- plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{x}\langle \partial^2_x\phi\rangle_y$'); hold on;
- plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{y}\langle \partial^2_x\phi\rangle_x$'); hold on;
- plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle \partial^2_x\phi\rangle_{x,y}$'); hold on;
- plot(Ts3D(its2D:ite2D),ones(ite2D-its2D+1,1)*shear_infty_avg, '-k',...
- 'DisplayName',['$s^{\infty} = $',num2str(shear_infty_avg),'$\pm$',num2str(shear_infty_std)]);
- ylim([0,shear_infty_avg*5.0]); xlim([Ts0D(1),Ts0D(end)]);
- grid on; ylabel('Shear amp.');set(gca,'xticklabel',[]);% legend('show');
+ it0 = 1; itend = Ns3D;
+ trange = it0:itend;
+ pltx = @(x) x;%-x(1);
+ plty = @(x) x./max((x(1:end)));
+% plty = @(x) x(500:end)./max(squeeze(x(500:end)));
+ yyaxis left
+ plot(pltx(Ts3D(trange)),plty(Ephi_Z(trange)),'DisplayName',['Zonal, ',CONAME]); hold on;
+ ylim([1e-1, 1.5]); ylabel('$E_{Z}$')
+ yyaxis right
+ plot(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt0(trange)),'DisplayName','$k_p>0$');
+% semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt1(trange)),'DisplayName','$k_p>1$');
+% semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt2(trange)),'DisplayName','$k_p>2$');
+ % semilogy(pltx(Ts0D),plty(PGAMMA_RI),'DisplayName',['$\Gamma_x$, ',CONAME]);
+% legend('Location','southeast')
+ xlim([Ts3D(it0), Ts3D(itend)]);
+ ylim([1e-6, 1.0]); ylabel('$E_{NZ}$')
+ xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
subplot(313)
[TY,TX] = meshgrid(x,Ts3D);
- pclr = pcolor(TX,TY,squeeze((mean(dx2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_x^2\phi\rangle_y$') %colorbar;
- caxis([-1 1]);
+ pclr = pcolor(TX,TY,squeeze((mean(dx2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_x^2\phi\rangle_y$') %colorbar;
+ caxis([-3 3]);
xlabel('$t c_s/R$'), ylabel('$x/\rho_s$'); colormap(bluewhitered(256))
-save_figure
\ No newline at end of file
+save_figure
diff --git a/matlab/plots/plot_space_time_diagrams.m b/matlab/plots/plot_space_time_diagrams.m
index 4fb784998788f35fb1aa13ef5666851f51cd9477..f82be30b4dedc722b412dd88ddba3281f0b2c14a 100644
--- a/matlab/plots/plot_space_time_diagrams.m
+++ b/matlab/plots/plot_space_time_diagrams.m
@@ -13,7 +13,7 @@ fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position', [100, 10
caxis([0.0,cmax]); c = colorbar; c.Label.String ='\langle\Gamma_{x}\rangle_{z}';
xticks([]); ylabel('$x/\rho_s$')
% legend('Radial part. transport $\langle n_i\partial_z\phi\rangle_z$')
- title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
+ title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\kappa_N=$',num2str(K_N),...
', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
' $\mu_{hd}=$',num2str(MU)]);
subplot(212)
@@ -23,5 +23,5 @@ fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position', [100, 10
caxis([-fieldmax,fieldmax]); c = colorbar; c.Label.String ='\langle v_{E\times B,z}\rangle_z';
colormap(bluewhitered)
xlabel('$t c_s/R$'), ylabel('$x/\rho_s$')
-% legend('Zonal flow $\langle \partial_r\phi\rangle_z$')
-save_figure
\ No newline at end of file
+% legend('Zonal flow $\langle \partial_r\phi\rangle_z$')
+save_figure
diff --git a/matlab/plots/plot_time_evolution_and_gr.m b/matlab/plots/plot_time_evolution_and_gr.m
index 17dbbb666db7a7d84171f9019ada51459bb86c64..b1f038b3af7392d4cf0b4c7fab62f3b9a741d2be 100644
--- a/matlab/plots/plot_time_evolution_and_gr.m
+++ b/matlab/plots/plot_time_evolution_and_gr.m
@@ -1,10 +1,10 @@
fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM),'_',PARAMS];
set(gcf, 'Position', [100, 100, 900, 800])
-subplot(111);
- suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
+subplot(111);
+ suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\kappa_N=$',num2str(K_N),...
', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
' $\mu_{hd}=$',num2str(MU)]);
- subplot(421);
+ subplot(421);
for ip = 1:Pe_max
for ij = 1:Je_max
plt = @(x) squeeze(x(ip,ij,:));
@@ -45,7 +45,7 @@ subplot(111);
ylim([0,max(max(g_I(1,:,:)))]); xlim([0,max(ky)]);
shearplot = 426; phiplot = 428;
else
- shearplot = 223; phiplot = 224;
+ shearplot = 223; phiplot = 224;
end
subplot(shearplot)
plt = @(x) mean(x,1);
@@ -55,7 +55,7 @@ subplot(111);
plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{r}\langle s \rangle_z$'); hold on;
plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{z}\langle s \rangle_r$'); hold on;
plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle s \rangle_{r,z}$'); hold on;
- grid on; xlabel('$t c_s/R$'); ylabel('$shear$');
+ grid on; xlabel('$t c_s/R$'); ylabel('$shear$');
subplot(phiplot)
clr = line_colors(min(ip,numel(line_colors(:,1))),:);
lstyle = line_styles(min(ij,numel(line_styles)));
@@ -64,4 +64,4 @@ subplot(111);
plot(Ts3D,plt(phi_avgx_maxy),'DisplayName','$\max_{z}\langle\phi\rangle_r$'); hold on;
plot(Ts3D,plt(phi_avgx_avgy),'DisplayName','$\langle\phi\rangle_{r,z}$'); hold on;
grid on; xlabel('$t c_s/R$'); ylabel('$E.S. pot$');
-save_figure
\ No newline at end of file
+save_figure
diff --git a/matlab/post_processing.m b/matlab/post_processing.m
index 1a74cb29e04777d1391abd142da91d5887f90484..14e8b6736d6d7c0264ec5d49173964997a9523a3 100644
--- a/matlab/post_processing.m
+++ b/matlab/post_processing.m
@@ -18,9 +18,9 @@ KPERP2 = KY.^2+KX.^2;
[KZ_XZ,KX_XZ] = meshgrid(z,kx);
[KZ_YZ,KY_YZ] = meshgrid(z,ky);
-Lk = max(Lkx,Lky);
-Nx = max(Nkx,Nky); Ny = Nx; Nz = numel(z);
-dx = 2*pi/Lk; dy = 2*pi/Lk; dz = 2*pi/Nz;
+Nx = 2*(Nkx-1); Ny = Nky; Nz = numel(z);
+Lx = 2*pi/dkx; Ly = 2*pi/dky;
+dx = Lx/Nx; dy = Ly/Ny; dz = 2*pi/Nz;
x = dx*(-Nx/2:(Nx/2-1)); Lx = max(x)-min(x);
y = dy*(-Ny/2:(Ny/2-1)); Ly = max(y)-min(y);
z = dz * (1:Nz);
@@ -77,7 +77,7 @@ for it = 1:numel(Ts3D)
temp_i (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_),Nx,Ny)));
Z_T_e (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_.*(KY==0)),Nx,Ny)));
Z_T_i (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_.*(KY==0)),Nx,Ny)));
- end
+ end
end
end
@@ -125,15 +125,15 @@ for it = 1:numel(Ts3D) % Loop over 2D aX_XYays
phi_avgx_maxy(iz,it) = max(mean(squeeze(phi(:,:,iz,it))));
phi_maxx_avgy(iz,it) = mean( max(squeeze(phi(:,:,iz,it))));
phi_avgx_avgy(iz,it) = mean(mean(squeeze(phi(:,:,iz,it))));
-
+
if(W_DENS)
Gamma_x(:,:,iz,it) = -dens_i(:,:,iz,it).*dyphi(:,:,iz,it);
Gamma_y(:,:,iz,it) = dens_i(:,:,iz,it).*dxphi(:,:,iz,it);
end
-
+
if(W_TEMP)
Q_x(:,:,iz,it) = -temp_e(:,:,iz,it).*dyphi(:,:,iz,it);
- Q_y(:,:,iz,it) = temp_i(:,:,iz,it).*dxphi(:,:,iz,it);
+ Q_y(:,:,iz,it) = temp_i(:,:,iz,it).*dxphi(:,:,iz,it);
end
shear_maxx_maxy(iz,it) = max( max(squeeze(-(dx2phi(:,:,iz,it)))));
@@ -171,7 +171,7 @@ for ikx = 1:Nkx
end
[gmax_I,ikmax_I] = max(max(g_I(1,:,:),[],2),[],3);
kmax_I = abs(ky(ikmax_I));
-Bohm_transport = ETAB/ETAN*gmax_I/kmax_I^2;
+Bohm_transport = K_N*gmax_I/kmax_I^2;
%% Compute secondary instability growth rate
disp('- growth rate')
@@ -192,3 +192,25 @@ for ikx = 1:Nkx
end
[gmax_II,ikmax_II] = max(max(g_II(1,:,:),[],2),[],3);
kmax_II = abs(kx(ikmax_II));
+
+%% zonal vs nonzonal energies for phi(t)
+Ephi_Z = zeros(1,Ns3D);
+Ephi_NZ_kgt0 = zeros(1,Ns3D);
+Ephi_NZ_kgt1 = zeros(1,Ns3D);
+Ephi_NZ_kgt2 = zeros(1,Ns3D);
+high_k_phi = zeros(1,Ns3D);
+for it = 1:numel(Ts3D)
+% Ephi_NZ(it) = sum(sum(((KY~=0).*abs(PHI(:,:,1,it)).^2)));
+% Ephi_Z(it) = sum(sum(((KY==0).*abs(PHI(:,:,1,it)).^2)));
+ [amp,ikzf] = max(abs((kx~=0).*PHI(:,1,1,it)));
+% Ephi_NZ(it) = sum(sum(((KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)));
+ Ephi_NZ_kgt0(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>0.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2))));
+ Ephi_NZ_kgt1(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>1.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2))));
+ Ephi_NZ_kgt2(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>2.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2))));
+% Ephi_Z(it) = kx(ikzf)^2*abs(PHI(ikzf,1,1,it)).^2;
+ Ephi_Z(it) = squeeze(sum(sum(((KX~=0).*(KY==0).*(KX.^2).*abs(PHI(:,:,1,it)).^2))));
+% Ephi_NZ(it) = sum(sum(((KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)))-Ephi_Z(it);
+ high_k_phi(it) = squeeze(abs(PHI(18,18,1,it)).^2); % kperp sqrt(2)
+% high_k_phi(it) = abs(PHI(40,40,1,it)).^2;% kperp 3.5
+
+end
diff --git a/matlab/setup.m b/matlab/setup.m
index 2587142b70d8088d9d2ed4c290e4981f8f6c61d8..e272a4f33cd119fc90ebf712e7033925291f8434 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -5,10 +5,10 @@ GRID.pmaxe = PMAXE; % Electron Hermite moments
GRID.jmaxe = JMAXE; % Electron Laguerre moments
GRID.pmaxi = PMAXI; % Ion Hermite moments
GRID.jmaxi = JMAXI; % Ion Laguerre moments
-GRID.Nx = N; % x grid resolution
-GRID.Lx = L; % x length
-GRID.Ny = N * (1-KXEQ0) + KXEQ0; % y ''
-GRID.Ly = L * (1-KXEQ0); % y ''
+GRID.Nx = Nx; % x grid resolution
+GRID.Lx = Lx; % x length
+GRID.Ny = Ny; % y ''
+GRID.Ly = Ly; % y ''
GRID.Nz = Nz; % z resolution
GRID.q0 = q0; % q factor
GRID.shear = shear; % shear
@@ -27,16 +27,18 @@ MODEL.nu = NU; % hyper diffusive coefficient nu for HW
MODEL.tau_e = TAU;
MODEL.tau_i = TAU;
% mass ratio sqrt(m_a/m_i)
-MODEL.sigma_e = 0.0233380;
+MODEL.sigma_e = SIGMA_E;
MODEL.sigma_i = 1.0;
% charge q_a/e
MODEL.q_e =-1.0;
MODEL.q_i = 1.0;
if MODEL.q_e == 0; SIMID = [SIMID,'_i']; end;
% gradients L_perp/L_x
-MODEL.eta_n = ETAN; % source term kappa for HW
-MODEL.eta_T = ETAT; % Temperature
-MODEL.eta_B = ETAB; % Magnetic
+MODEL.K_n = K_N; % source term kappa for HW
+MODEL.K_T = K_T; % Temperature
+MODEL.K_E = K_E; % Electric
+MODEL.GradB = GRADB; % Magnetic gradient
+MODEL.CurvB = CURVB; % Magnetic curvature
MODEL.lambdaD = LAMBDAD;
% if A0KH ~= 0; SIMID = [SIMID,'_Nz_',num2str(L/2/pi*KX0KH),'_A_',num2str(A0KH)]; end;
% Time integration and intialization parameters
@@ -79,27 +81,45 @@ if (CLOS == 0); CLOSNAME = 'Trunc.';
elseif(CLOS == 1); CLOSNAME = 'Clos. 1';
elseif(CLOS == 2); CLOSNAME = 'Clos. 2';
end
+% Hermite-Laguerre degrees naming
if (PMAXE == PMAXI) && (JMAXE == JMAXI)
- degngrad = ['P_',num2str(PMAXE),'_J_',num2str(JMAXE)];
+ HLdeg_ = ['_',num2str(PMAXE+1),'x',num2str(JMAXE+1)];
else
- degngrad = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),...
- '_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI)];
+ HLdeg_ = ['_Pe_',num2str(PMAXE+1),'_Je_',num2str(JMAXE+1),...
+ '_Pi_',num2str(PMAXI+1),'_Ji_',num2str(JMAXI+1)];
end
-degngrad = [degngrad,'_eta_',num2str(ETAB/ETAN),'_nu_%0.0e_',...
- CONAME,'_mu_%0.0e'];
-
-degngrad = sprintf(degngrad,[NU,MU]);
-if ~NON_LIN; degngrad = ['lin_',degngrad]; end
-if (Nz == 1)
- resolution = [num2str(GRID.Nx),'x',num2str(GRID.Ny/2),'_'];
- gridname = ['L_',num2str(L),'_'];
+% temp. dens. drives
+drives_ = [];
+if abs(K_N) > 0; drives_ = [drives_,'_kN_',num2str(K_N)]; end;
+if abs(K_T) > 0; drives_ = [drives_,'_kT_',num2str(K_T)]; end;
+% collision
+coll_ = ['_nu_%0.0e_',CONAME];
+coll_ = sprintf(coll_,NU);
+% nonlinear
+lin_ = [];
+if ~NON_LIN; lin_ = '_lin'; end
+% resolution and boxsize
+res_ = [num2str(GRID.Nx),'x',num2str(GRID.Ny)];
+if (Lx ~= Ly)
+ geo_ = ['_Lx_',num2str(Lx),'_Ly_',num2str(Ly)];
else
- resolution = [num2str(GRID.Nx),'x',num2str(GRID.Ny/2),'x',num2str(GRID.Nz),'_'];
- gridname = ['L_',num2str(L),'_q0_',num2str(q0),'_'];
+ geo_ = ['_L_',num2str(Lx)];
+end
+if (Nz > 1) %3D case
+ res_ = [res_,'x',num2str(GRID.Nz)];
+ if abs(q0) > 0
+ geo_ = [geo_,'_q0_',num2str(q0)];
+ end
+ if abs(eps) > 0
+ geo_ = [geo_,'_e_',num2str(eps)];
+ end
+ if abs(shear) > 0
+ geo_ = [geo_,'_s_',num2str(shear)];
+ end
end
-if (exist('PREFIX','var') == 0); PREFIX = []; end;
-if (exist('SUFFIX','var') == 0); SUFFIX = []; end;
-PARAMS = [PREFIX,resolution,gridname,degngrad,SUFFIX];
+% put everything together in the param character chain
+u_ = '_'; % underscore variable
+PARAMS = [res_,HLdeg_,geo_,drives_,coll_,lin_];
BASIC.RESDIR = [SIMDIR,PARAMS,'/'];
BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',SIMDIR(4:end),PARAMS,'/'];
BASIC.PARAMS = PARAMS;
@@ -116,7 +136,6 @@ OUTPUTS.nsave_1d = -1;
OUTPUTS.nsave_2d = floor(1.0/SPS2D/DT);
OUTPUTS.nsave_3d = floor(1.0/SPS3D/DT);
OUTPUTS.nsave_5d = floor(1.0/SPS5D/DT);
-OUTPUTS.nsave_cp = floor(1.0/SPSCP/DT);
if W_DOUBLE; OUTPUTS.write_doubleprecision = '.true.'; else; OUTPUTS.write_doubleprecision = '.false.';end;
if W_GAMMA; OUTPUTS.write_gamma = '.true.'; else; OUTPUTS.write_gamma = '.false.';end;
if W_HF; OUTPUTS.write_hf = '.true.'; else; OUTPUTS.write_hf = '.false.';end;
@@ -144,7 +163,7 @@ MAKE = 'cd ..; make; cd wk';
% system(MAKE);
%%
disp(['Set up ',SIMID]);
-disp([resolution,gridname,degngrad]);
+disp([res_,geo_,HLdeg_]);
if JOB2LOAD>=0
disp(['- restarting from JOBNUM = ',num2str(JOB2LOAD)]); else
disp(['- starting from T = 0']);
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 19b918c376ac2ca883b08cc613ca58e06814af43..57873487a7619ca979d19fb1059acbbb446f92cb 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -59,9 +59,11 @@ fprintf(fid,[' sigma_e = ', num2str(MODEL.sigma_e),'\n']);
fprintf(fid,[' sigma_i = ', num2str(MODEL.sigma_i),'\n']);
fprintf(fid,[' q_e = ', num2str(MODEL.q_e),'\n']);
fprintf(fid,[' q_i = ', num2str(MODEL.q_i),'\n']);
-fprintf(fid,[' eta_n = ', num2str(MODEL.eta_n),'\n']);
-fprintf(fid,[' eta_T = ', num2str(MODEL.eta_T),'\n']);
-fprintf(fid,[' eta_B = ', num2str(MODEL.eta_B),'\n']);
+fprintf(fid,[' K_n = ', num2str(MODEL.K_n),'\n']);
+fprintf(fid,[' K_T = ', num2str(MODEL.K_T),'\n']);
+fprintf(fid,[' K_E = ', num2str(MODEL.K_E),'\n']);
+fprintf(fid,[' GradB = ', num2str(MODEL.GradB),'\n']);
+fprintf(fid,[' CurvB = ', num2str(MODEL.CurvB),'\n']);
fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']);
fprintf(fid,'/\n');
diff --git a/matlab/write_sbash_marconi.m b/matlab/write_sbash_marconi.m
index 481b72548d9ad8c175cf22fa3b78372606813185..28a55760a6c3c80fcf0165fa208a89c21ab8bdf2 100644
--- a/matlab/write_sbash_marconi.m
+++ b/matlab/write_sbash_marconi.m
@@ -15,7 +15,7 @@ fprintf(fid,[...
'rm $HOME/HeLaZ/wk/fort*.90\n',...
...
SBATCH_CMD,...
-'echo tail -f $CINECA_SCRATCH/HeLaZ',BASIC.RESDIR(3:end),'out']);
+'echo tail -f $CINECA_SCRATCH/HeLaZ',BASIC.RESDIR(3:end)]);
fclose(fid);
system(['cp setup_and_run.sh ',BASIC.RESDIR,'/.']);
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index e661b401c309ccd1600d8f4f3a28a0a737f82172..5f32c1b52d389066964fb431353d3649a1030088 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -4,29 +4,24 @@ outfile ='';
%% Directory of the simulation
if 1% Local results
outfile ='';
-outfile ='';
-outfile ='';
-outfile ='';
-outfile ='';
-outfile ='';
-outfile ='';
-outfile ='';
-outfile ='test_even_p/100x50_L_100_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_0e+00';
+outfile ='simulation_A/1024x256_3x2_L_120_kN_1.6667_nu_1e-01_DGGK';
+% outfile ='Linear_Device/64x64x20_5x2_Lx_20_Ly_150_q0_25_kN_0.24_kT_0.03_nu_1e-02_DGGK';
BASIC.RESDIR = ['../results/',outfile,'/'];
BASIC.MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
system(['mkdir -p ',BASIC.MISCDIR]);
CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
system(CMD);
else% Marconi results
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/500x500_L_120_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_0e+00/out.txt';
BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
end
%% Load the results
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 20;
+JOBNUMMIN = 02; JOBNUMMAX = 02;
% JOBNUMMIN = 07; JOBNUMMAX = 20; % For CO damping sim A
compile_results %Compile the results from first output found to JOBNUMMAX if existing
@@ -45,7 +40,7 @@ end
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 = 1000; TAVG_1 = 2000; % Averaging times duration
+TAVG_0 = 1500; TAVG_1 = 4000; % Averaging times duration
plot_radial_transport_and_shear
end
@@ -59,7 +54,7 @@ end
if 0
%% |phi_k|^2 spectra (Kobayashi 2015 fig 3)
% tstart = 0.8*Ts3D(end); tend = Ts3D(end); % Time window
-tstart = 2500; tend = 2500;
+tstart = 1000; tend = 4000;
% tstart = 10000; tend = 12000;
% Chose the field to plot
% FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
@@ -74,20 +69,20 @@ if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
t0 =000; iz = 1; ix = 1; iy = 1;
-skip_ =1; DELAY = 2e-3*skip_;
+skip_ =1; FPS = 30; DELAY = 1/FPS;
[~, it03D] = min(abs(Ts3D-t0)); FRAMES_3D = it03D:skip_:numel(Ts3D);
[~, it05D] = min(abs(Ts5D-t0)); FRAMES_5D = it05D:skip_:numel(Ts5D);
T = Ts3D; FRAMES = FRAMES_3D;
-INTERP = 0; NORMALIZED = 0; CONST_CMAP = 0; BWR =1;% Gif options
+INTERP = 0; NORMALIZED = 1; CONST_CMAP = 0; BWR =1;% Gif options
% Field to plot
% FIELD = dens_i; NAME = 'ni'; FIELDNAME = 'n_i';
% FIELD = dens_i-Z_n_i; NAME = 'ni_NZ';FIELDNAME = 'n_i^{NZ}';
% FIELD = temp_i; NAME = 'Ti'; FIELDNAME = 'n_i';
% FIELD = temp_i-Z_T_i; NAME = 'Ti_NZ';FIELDNAME = 'T_i^{NZ}';
% FIELD = ne00; NAME = 'ne00'; FIELDNAME = 'n_e^{00}';
-% FIELD = ni00; NAME = 'ni00'; FIELDNAME = 'n_i^{00}';
-FIELD = phi; NAME = 'phi'; FIELDNAME = '\phi';
-% FIELD = Z_phi; NAME = 'Zphi'; FIELDNAME = '\phi_Z';
+FIELD = ni00; NAME = 'ni00'; FIELDNAME = 'n_i^{00}';
+% FIELD = phi; NAME = 'phi'; FIELDNAME = '\phi';
+% FIELD = phi-Z_phi; NAME = 'NZphi'; FIELDNAME = '\phi_{NZ}';
% FIELD = Gamma_x; NAME = 'Gamma_x'; FIELDNAME = '\Gamma_x';
% Sliced
@@ -95,8 +90,9 @@ FIELD = phi; NAME = 'phi'; FIELDNAME = '\phi';
% plt = @(x) real(x( :,iy, :,:)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z';
plt = @(x) real(x( :, :,iz,:)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y';
-% % K-space
-% % FIELD = PHI; NAME = 'PHI'; FIELDNAME = '\tilde \phi';
+% K-space
+% FIELD = PHI; NAME = 'PHI'; FIELDNAME = '\tilde \phi';
+% FIELD = Ne00; NAME = 'Ne00'; FIELDNAME = 'N_e^{00}';
% FIELD = Ni00; NAME = 'Ni00'; FIELDNAME = 'N_i^{00}';
% plt = @(x) fftshift((abs(x( :, :,1,:))),2); X = fftshift(KX,2); Y = fftshift(KY,2); XNAME = 'k_x'; YNAME = 'k_y';
@@ -111,8 +107,8 @@ FIELD = squeeze(plt(FIELD));
GIFNAME = [NAME,sprintf('_%.2d',JOBNUM),'_',PARAMS];
% Create movie (gif or mp4)
-create_gif
-% create_mov
+% create_gif
+create_mov
end
if 0
@@ -120,7 +116,7 @@ if 0
% Chose the field to plot
% FIELD = ni00; FNAME = 'ni00'; FIELDLTX = 'n_i^{00}';
-% FIELD = ne00; FNAME = 'ne00'; FIELDLTX = 'n_e^{00}'
+FIELD = ne00; FNAME = 'ne00'; FIELDLTX = 'n_e^{00}'
% FIELD = dens_i; FNAME = 'ni'; FIELDLTX = 'n_i';
% FIELD = dens_e; FNAME = 'ne'; FIELDLTX = 'n_e';
% FIELD = dens_e-Z_n_e; FNAME = 'ne_NZ'; FIELDLTX = 'n_e^{NZ}';
@@ -129,12 +125,12 @@ if 0
% FIELD = temp_e; FNAME = 'Te'; FIELDLTX = 'T_e';
% FIELD = phi; FNAME = 'phi'; FIELDLTX = '\phi';
% FIELD = Z_phi-phi; FNAME = 'phi_NZ'; FIELDLTX = '\phi^{NZ}';
-FIELD = Z_phi; FNAME = 'phi_Z'; FIELDLTX = '\phi^{Z}';
+% FIELD = Z_phi; FNAME = 'phi_Z'; FIELDLTX = '\phi^{Z}';
% FIELD = Gamma_x; FNAME = 'Gamma_x'; FIELDLTX = '\Gamma_x';
% FIELD = dens_e-Z_n_e-(Z_phi-phi); FNAME = 'Non_adiab_part'; FIELDLTX = 'n_e^{NZ}-\phi^{NZ}';
% Chose when to plot it
-tf = [9800 10000];
+tf = 200:200:1000;
% Sliced
ix = 1; iy = 1; iz = 1;
@@ -170,13 +166,13 @@ if 0
% Chose the field to plot
% FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
-% FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
-FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
+FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
+% FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
% FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
% FIELD_ = fft2(dens_e); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'FFT_Dens_e'; FIELDLTX = '\tilde n_e';
% Chose when to plot it
-tf = 1000:500:3000;
+tf = 200:400:1400;
% tf = 8000;
% Sliced
@@ -195,6 +191,7 @@ plt_2 = @(x) (fftshift(x,2));
DATA = plt_2(squeeze(plt(FIELD,it)));
pclr = pcolor(fftshift(X,2),fftshift(Y,2),DATA); set(pclr, 'edgecolor','none');pbaspect([0.5 1 1])
caxis([0 1]*5e3);
+% caxis([-1 1]*5);
xlabel(latexize(XNAME)); ylabel(latexize(YNAME));
if(i_ > 1); set(gca,'ytick',[]); end;
title(sprintf('$t c_s/R=%.0f$',Ts3D(it)));
@@ -216,7 +213,7 @@ end
if 0
%% Radial shear profile
-tf = 3000+[900:20:1100];
+tf = 2000+[900:20:1100];
ymax = 0;
figure
for i_ = 1:numel(tf)
@@ -231,27 +228,9 @@ end
if 1
%% zonal vs nonzonal energies for phi(t)
-trange = 1:Ns3D;
-Ephi_Z = zeros(1,Ns3D);
-Ephi_NZ_kgt0 = zeros(1,Ns3D);
-Ephi_NZ_kgt1 = zeros(1,Ns3D);
-Ephi_NZ_kgt2 = zeros(1,Ns3D);
-high_k_phi = zeros(1,Ns3D);
-for it = 1:numel(Ts3D)
-% Ephi_NZ(it) = sum(sum(((KY~=0).*abs(PHI(:,:,1,it)).^2)));
-% Ephi_Z(it) = sum(sum(((KY==0).*abs(PHI(:,:,1,it)).^2)));
- [amp,ikzf] = max(abs((kx~=0).*PHI(:,1,1,it)));
-% Ephi_NZ(it) = sum(sum(((KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)));
- Ephi_NZ_kgt0(it) = sum(sum(((sqrt(KX.^2+KY.^2)>0.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)));
- Ephi_NZ_kgt1(it) = sum(sum(((sqrt(KX.^2+KY.^2)>1.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)));
- Ephi_NZ_kgt2(it) = sum(sum(((sqrt(KX.^2+KY.^2)>2.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)));
-% Ephi_Z(it) = kx(ikzf)^2*abs(PHI(ikzf,1,1,it)).^2;
- Ephi_Z(it) = sum(sum(((KX~=0).*(KY==0).*(KX.^2).*abs(PHI(:,:,1,it)).^2)));
-% Ephi_NZ(it) = sum(sum(((KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)))-Ephi_Z(it);
- high_k_phi(it) = abs(PHI(18,18,1,it)).^2; % kperp sqrt(2)
-% high_k_phi(it) = abs(PHI(40,40,1,it)).^2;% kperp 3.5
-end
-pltx = @(x) x-x(1);
+it0 = 01; itend = Ns3D;
+trange = it0:itend;
+pltx = @(x) x;%-x(1);
plty = @(x) x./max(squeeze(x));
fig = figure; FIGNAME = ['ZF_turb_energy_vs_time_',PARAMS];
set(gcf, 'Position', [100, 100, 1400, 500])
@@ -263,21 +242,15 @@ subplot(121)
semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt1(trange)),'DisplayName',['NZ, $k_p>1$, ',CONAME]);
semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt2(trange)),'DisplayName',['NZ, $k_p>2$, ',CONAME]);
% semilogy(pltx(Ts0D),plty(PGAMMA_RI),'DisplayName',['$\Gamma_x$, ',CONAME]);
- title('Energy'); legend('show')
+ title('Energy'); legend('Location','southeast')
+ xlim([Ts3D(it0), Ts3D(itend)]);
+ ylim([1e-3, 1.5])
xlabel('$t c_s/R$'); grid on;% xlim([0 500]);
subplot(122)
plot(plty(Ephi_Z(trange)),plty(Ephi_NZ_kgt0(trange)));
title('Phase space'); legend(CONAME)
xlabel('$E_Z$'); ylabel('$E_{NZ}$'); grid on;% xlim([0 500]);
-%
-% subplot(133)
-% % semilogy(pltx(Ts0D),plty(abs(PGAMMA_RI)*SCALE));
-% for ik = [10 20 30]
-% semilogy(pltx(Ts3D(trange)),plty(abs(squeeze(PHI(ik,ik,1,trange))).^2),'DisplayName',[CONAME,', kp=',num2str(sqrt(kx(ik)^2+ky(ik)^2))]); hold on
-% end
-% title('High kperp damping'); legend('show');
-% xlabel('$t c_s/R$'); grid on;% xlim([0 500]);
end
if 0
@@ -305,4 +278,35 @@ subplot(212);
% plot(Ts3D,squeeze(mflux_y_i),'DisplayName','Flux i y');
% plot(Ts3D,squeeze(mflux_y_o),'DisplayName','Flux o y');
legend('show'); grid on; xlim([Ts3D(1) Ts3D(end)]); %ylim([-0.1, 2]*mean(mflux_x_i))
+end
+
+if 0
+%% Zonal profiles (ky=0)
+
+% Chose the field to plot
+FIELD = Ne00.*conj(Ne00); FNAME = 'Ne002'; FIELDLTX = '|N_e^{00}|^2'
+% FIELD = Ni00.*conj(Ni00); FNAME = 'Ni002'; FIELDLTX = '|N_i^{00}|^2'
+% FIELD = abs(PHI); FNAME = 'absPHI'; FIELDLTX = '|\tilde\phi|';
+% FIELD = PHI.*conj(PHI); FNAME = 'PHI2'; FIELDLTX = '|\tilde\phi|^2';
+% FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
+% FIELD_ = fft2(dens_e); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'FFT_Dens_e'; FIELDLTX = '\tilde n_e';
+
+% Chose when to plot it
+tf = 1500:200:2500;
+% tf = 8000;
+
+% Sliced
+plt = @(x,it) x( 2:end, 1,1,it)./max(max(x( 2:end, 1,1,:))); XNAME = 'k_x';
+%
+TNAME = [];
+fig = figure; FIGNAME = ['Zonal_',FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 600,400])
+plt_2 = @(x) (fftshift(x,2));
+ for i_ = 1:numel(tf)
+ [~,it] = min(abs(Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(Ts3D(it))];
+ DATA = plt_2(squeeze(plt(FIELD,it)));
+ semilogy(kx(2:end),DATA,'-','DisplayName',sprintf('$t c_s/R=%.0f$',Ts3D(it))); hold on; grid on;
+ xlabel(latexize(XNAME));
+ end
+title(['$',FIELDLTX,'$ Zonal Spectrum']); legend('show');
+save_figure
end
\ No newline at end of file
diff --git a/wk/continue_multiple_runs_marconi.m b/wk/continue_multiple_runs_marconi.m
index ea1fa75c00f0a242920a366cc19587bda10be206..21cf32acecd84b28fe906c84c10119fc6a10f6ad 100644
--- a/wk/continue_multiple_runs_marconi.m
+++ b/wk/continue_multiple_runs_marconi.m
@@ -1,5 +1,5 @@
%% Paste the list of continue_run calls
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt')
+% continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt')
continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt')
%% Functions to modify preexisting fort.90 input file and launch on marconi
@@ -39,18 +39,10 @@ function [] = continue_run(outfilename)
fclose(fid);
% Find previous job2load
- if( numel(A{5}) == numel(' RESTART = .false.') )
- A{5} = sprintf(' RESTART = .true.');
- J2L = 0;
- else
- line = A{39};
- line = line(end-2:end);
- if(line(1) == '='); line = line(end); end;
- J2L = str2num(line)+1; JOB2LOAD=J2L;
- end
- % Change job 2 load in fort.90
- A{39} = [' job2load = ',num2str(J2L)];
- disp(A{39})
+ line = A{38};
+ J2L_old = str2num(line(17:end););
+ A{38} = [' job2load = ',num2str(J2L_old+1)];
+ disp(A{38})
% Change time step
line_= A{3};
dt_old = str2num(line_(13:end));
@@ -58,28 +50,28 @@ function [] = continue_run(outfilename)
% Increase endtime
A{4} = [' tmax = 20000'];
% Change collision operator
- line_= A{43};
+ line_= A{42};
CO_old = str2num(line_(13:end));
- A{43} = [' CO = ',num2str(2)];
+ A{42} = [' CO = ',num2str(2)];
% Put non linear term back
- A{45} = [' NL_CLOS = -1'];
+ A{44} = [' NL_CLOS = -1'];
% change HD
- line_= A{47};
+ line_= A{46};
mu_old = str2num(line_(13:end));
- A{47} = [' mu = ',num2str(0)];
+ A{46} = [' mu = ',num2str(0)];
% change L
- line_= A{14};
+ line_= A{13};
L_old = str2num(line_(12:end));
- A{14} = [' Lx = ',num2str(L_old)];
- A{16} = [' Ly = ',num2str(L_old)];
+ A{13} = [' Lx = ',num2str(L_old)];
+ A{15} = [' Ly = ',num2str(L_old)];
% change eta N
- line_= A{57};
+ line_= A{56};
etan_old = str2num(line_(13:end));
- A{57} = [' eta_n = ',num2str(etan_old)];
+ A{56} = [' eta_n = ',num2str(etan_old)];
% change eta B
- line_= A{59};
+ line_= A{58};
etab_old = str2num(line_(13:end));
- A{59} = [' eta_B = ',num2str(etab_old)];
+ A{58} = [' eta_B = ',num2str(etab_old)];
% Rewrite fort.90
fid = fopen('fort.90', 'w');
for i = 1:numel(A)
diff --git a/wk/kobayashi_2015_linear_results.m b/wk/kobayashi_2015_linear_results.m
new file mode 100644
index 0000000000000000000000000000000000000000..82f339125f78966b0d059746be9b56d8eaa292d0
--- /dev/null
+++ b/wk/kobayashi_2015_linear_results.m
@@ -0,0 +1,19 @@
+% results for blue triangles in Kobayashi 2015
+% x_ = [0.06 0.4 0.6];
+% y_ = [0.0 0.1 0.0];
+% plot(x_,y_,'^b','DisplayName','Kobayashi et al. 2015');
+
+
+% results for red squares in Kobayashi 2015
+% x_ = [0.065 0.15 0.35];
+% y_ = [0.0 0.15 0.0];
+% plot(x_,y_,'sr','DisplayName','Kobayashi et al. 2015');
+
+
+% results for green triangles in Kobayashi 2015
+x_ = [0.06 0.5 1.8];
+y_ = [0.0 0.2 0.0];
+plot(x_,y_,'dg','DisplayName','Kobayashi et al. 2015');
+
+
+
diff --git a/wk/linear_study.m b/wk/linear_1D_entropy_mode.m
similarity index 58%
rename from wk/linear_study.m
rename to wk/linear_1D_entropy_mode.m
index eed525c5b6e8ed10e735cf5aafa17f3598a32ddd..0408124767362a246a30463007ac544a6301a1bc 100644
--- a/wk/linear_study.m
+++ b/wk/linear_1D_entropy_mode.m
@@ -1,5 +1,4 @@
-for CO = [3]
- RUN = 1; % To run or just to load
+RUN = 1; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -7,43 +6,39 @@ default_plots_options
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.1; % Collision frequency
+NU = 0.01; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-ETAB = 1.0;
-ETAN = 1./0.6; % Density gradient
-ETAT = 0.0; % Temperature gradient
-NU_HYP = 0.0; % Hyperdiffusivity coefficient
-LAMBDAD = 0.0;
-NOISE0 = 1.0e-5; % Init noise amplitude
-BCKGD0 = 0.0; % Init background
+K_N = 1/0.6; % Density gradient drive
+K_T = 0.0; % Temperature '''
+K_E = 0.0; % Electrostat '''
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
%% GRID PARAMETERS
-N = 150; % Frequency gridpoints (Nkx = N/2)
-L = 100; % Size of the squared frequency domain
-KXEQ0 = 1; % put kx = 0
-MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
-MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
+Nx = 82; % real space x-gridpoints
+Ny = 1; % '' y-gridpoints
+Lx = 150; % Size of the squared frequency domain
+Ly = 0; % Size of the squared frequency domain
+Nz = 1; % number of perpendicular planes (parallel grid)
+q0 = 1.0; % safety factor
+shear = 0.0; % magnetic shear
+eps = 0.0; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 200; % Maximal time unit
+TMAX = 300; % 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
-SPS5D = 2; % Sampling per time unit for 5D arrays
+SPS3D = 1; % Sampling per time unit for 2D arrays
+SPS5D = 1; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
-% SIMID = 'v3.6_kobayashi_lin'; % Name of the simulation
-% SIMID = 'v3.2_CO_damping'; % Name of the simulation
-% SIMID = 'CO_Patchwork_damping'; % Name of the simulation
-SIMID = 'test_even_p'; % Name of the simulation
-% SIMID = 'v3.2_entropy_mode_linear'; % Name of the simulation
-NON_LIN = 0 *(1-KXEQ0); % activate non-linearity (is cancelled if KXEQ0 = 1)
+SIMID = 'test_4.1'; % Name of the simulation
+NON_LIN = 0; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
-% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle, 4 : Full Couloumb ; +/- for GK/DK)
-% CO = 1;
+% (0:L.Bernstein, 1:Dougherty, 2:Sugama, 3:Pitch angle, 4:Full Couloumb ; +/- for GK/DK)
+CO = 3;
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))
-NL_CLOS = 0; % nonlinear closure model (0: =0 nmax = jmax, 1: nmax = jmax-j, >1 : nmax = NL_CLOS)
+NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
INIT_PHI= 0; % Start simulation with a noisy phi
%% OUTPUTS
@@ -56,32 +51,34 @@ W_NAPJ = 1; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
-kmax = N*pi/L;% Highest fourier mode
+kmax = Nx*pi/Lx;% Highest fourier mode
+NU_HYP = 0.0; % Hyperdiffusivity coefficient
MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
-Nz = 1; % number of perpendicular planes (parallel grid)
-q0 = 1.0; % safety factor
-shear = 0.0; % magnetic shear
-eps = 0.0; % inverse aspect ratio
INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0;
+MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
+MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
+LAMBDAD = 0.0;
+NOISE0 = 1.0e-5; % Init noise amplitude
+BCKGD0 = 0.0; % Init background
+GRADB = 1.0;
+CURVB = 1.0;
%% PARAMETER SCANS
if 1
%% Parameter scan over PJ
% PA = [2 4];
% JA = [1 2];
-PA = [8];
-JA = [4];
+PA = [4];
+JA = [2];
DTA= DT*ones(size(JA));%./sqrt(JA);
% DTA= DT;
mup_ = MU_P;
muj_ = MU_J;
Nparam = numel(PA);
param_name = 'PJ';
-gamma_Ni00 = zeros(Nparam,floor(N/2)+1);
-gamma_Nipj = zeros(Nparam,floor(N/2)+1);
-gamma_phi = zeros(Nparam,floor(N/2)+1);
-% Ni00_ST = zeros(Nparam,floor(N/2)+1,TMAX/SPS3D);
-% PHI_ST = zeros(Nparam,floor(N/2)+1,TMAX/SPS3D);
+gamma_Ni00 = zeros(Nparam,floor(Nx/2)+1);
+gamma_Nipj = zeros(Nparam,floor(Nx/2)+1);
+gamma_phi = zeros(Nparam,floor(Nx/2)+1);
for i = 1:Nparam
% Change scan parameter
PMAXE = PA(i); PMAXI = PA(i);
@@ -99,28 +96,28 @@ for i = 1:Nparam
%%
filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5'];
load_results
- for ikx = 1:N/2+1
- %find if there is a tail of 0 (highest damped mode)
-
+ for ikx = 1:Nx/2+1
tend = max(Ts3D(abs(Ni00(ikx,1,1,:))~=0));
- tstart = 0.8*tend;
+ tstart = 0.6*tend;
[~,itstart] = min(abs(Ts3D-tstart));
[~,itend] = min(abs(Ts3D-tend));
- gamma_Ni00(i,ikx) = (LinearFit_s(Ts3D(itstart:itend)',(squeeze(abs(Ni00(ikx,1,1,itstart:itend))))));
- gamma_phi (i,ikx) = (LinearFit_s(Ts3D(itstart:itend)',(squeeze(abs(PHI (ikx,1,1,itstart:itend))))));
-% Ni00_ST(i,ikx,:) = squeeze((Ni00(ikx,1,1,1:TMAX/SPS3D)));
-% PHI_ST(i,ikx,:) = squeeze((PHI (ikx,1,1,1:TMAX/SPS3D)));
- end
- tend = Ts5D(end); tstart = 0.4*tend;
- [~,itstart] = min(abs(Ts5D-tstart));
- [~,itend] = min(abs(Ts5D-tend));
- for ikx = 1:N/2+1
- gamma_Nipj(i,ikx) = LinearFit_s(Ts5D(itstart:itend)',squeeze(max(max(abs(Nipj(:,:,ikx,1,1,itstart:itend)),[],1),[],2)));
+ trange = itstart:itend;
+ % exp fit on moment 00
+ X_ = Ts3D(trange); Y_ = squeeze(abs(Ni00(ikx,1,1,trange)));
+ gamma_Ni00(i,ikx) = LinearFit_s(X_,Y_);
+ % exp fit on phi
+ X_ = Ts3D(trange); Y_ = squeeze(abs(PHI(ikx,1,1,trange)));
+ gamma_phi (i,ikx) = LinearFit_s(X_,Y_);
end
gamma_Ni00(i,:) = real(gamma_Ni00(i,:));% .* (gamma_Ni00(i,:)>=0.0));
gamma_Nipj(i,:) = real(gamma_Nipj(i,:));% .* (gamma_Nipj(i,:)>=0.0));
- % Clean output
-% system(['rm -r ',BASIC.RESDIR]);
+ if 1
+ %% Fit verification
+ figure;
+ for i = 1:4:Nx/2+1
+ X_ = Ts3D(:); Y_ = squeeze(abs(Ni00(i,1,1,:)));
+ semilogy(X_,Y_,'DisplayName',['k=',num2str(kx(i))]); hold on;
+ end
end
if 1
@@ -132,15 +129,15 @@ plt = @(x) x;
for i = 1:Nparam
clr = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:);
linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1);
- plot(plt(SCALE*kx(1:numel(kx))),plt(gamma_Ni00(i,1:end)),...
+ plot(plt(SCALE*kx),plt(gamma_phi(i,1:end)),...
'Color',clr,...
'LineStyle',linestyle{1},'Marker','^',...
-...% 'DisplayName',['$\eta=',num2str(ETAB/ETAN),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, $P=',num2str(PA(i)),'$, $J=',num2str(JA(i)),'$']);
+...% 'DisplayName',['$\kappa_N=',num2str(K_N),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, $P=',num2str(PA(i)),'$, $J=',num2str(JA(i)),'$']);
'DisplayName',[CONAME,', $P,J=',num2str(PA(i)),',',num2str(JA(i)),'$']);
hold on;
end
- grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(N_i^{00})L_\perp/c_s$'); xlim([0.0,max(kx)]);
- title(['$\eta=',num2str(ETAB/ETAN),'$, $\nu_{',CONAME,'}=',num2str(NU),'$'])
+ grid on; xlabel('$k_y\rho_s^{R}$'); ylabel('$\gamma(\phi)L_\perp/c_s$'); xlim([0.0,max(kx)]);
+ title(['$\kappa_N=',num2str(K_N),'$, $\nu_{',CONAME,'}=',num2str(NU),'$'])
% title(['$\nabla N = 0$', ', $\nu=',num2str(NU),'$'])
legend('show'); %xlim([0.01,10])
saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']);
@@ -155,11 +152,4 @@ if 0
% pclr = pcolor(XT,YT,squeeze(abs(Ni00_ST(1,:,:)))); set(pclr, 'edgecolor','none'); colorbar;
semilogy(Ts3D(1:TMAX/SPS3D),squeeze(abs(PHI_ST(1,50:5:100,:))));
end
-if 0
-%% Trajectories of some modes
-figure;
-for i = 1:10:N/2+1
- semilogy(Ts3D,squeeze(abs(Ne00(i,1,1,:))),'DisplayName',['k=',num2str(kx(i))]); hold on;
-end
-end
-end
+end
\ No newline at end of file
diff --git a/wk/load_multiple_outputs_marconi.m b/wk/load_multiple_outputs_marconi.m
index 0738cad9421e410b82e9d123eb6b13a1d2869ec2..1fdd3b8608236ae4a056034abf7852fc50a814a0 100644
--- a/wk/load_multiple_outputs_marconi.m
+++ b/wk/load_multiple_outputs_marconi.m
@@ -1,4 +1,6 @@
addpath(genpath('../matlab')) % ... add
-%% Paste the list of simulation results to load
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/500x500_L_120_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_0e+00/')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/500x500_L_120_P_2_J_1_eta_0.6_nu_1e-01_DGGK_mu_0e+00/')
diff --git a/wk/local_run.m b/wk/local_run.m
index 3a3472529a42e1555a946eadd641983c3c36a8e1..0120ca24cf0adabcc1671a10acc075916c2d65bf 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -5,24 +5,31 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
-ETAN = 1.0/0.6; % Density gradient drive (R/Ln)
+K_N = 1/0.6; % Density gradient drive (R/Ln)
+K_T = 0.00; % Temperature gradient
+K_E = 0.00; % Electrostat gradient
NU_HYP = 0.0;
-%% GRID AND GEOMETRY PARAMETERS
-N = 100; % Frequency gridpoints (Nkx = N/2)
-L = 100; % Size of the squared frequency domain
-Nz = 1; % number of perpendicular planes (parallel grid)
-q0 = 1.0; % safety factor
-shear = 0.0; % magnetic shear
-eps = 0.0; % inverse aspect ratio
-P = 4;
-J = 2;
+%% GRID PARAMETERS
+Nx = 1024; % Spatial radial resolution ( = 2x radial modes)
+Lx = 120; % Radial window size
+Ny = 256; % Spatial azimuthal resolution (= azim modes)
+Ly = 120; % Azimuthal window size
+Nz = 1; % number of perpendicular planes (parallel grid)
+q0 = 1.0; % safety factor (Lz = 2*pi*q0)
+P = 2;
+J = 1;
+%% GEOMETRY PARAMETERS
+shear = 0.0; % magnetic shear
+eps = 0.0; % inverse aspect ratio (controls parallel magnetic gradient)
+gradB = 0.0; % Magnetic gradient
+curvB = 0.0; % Magnetic curvature
%% TIME PARAMETERS
-TMAX = 100; % Maximal time unit
+TMAX = 90; % Maximal time unit
DT = 2e-2; % Time step
SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS3D = 1; % Sampling per time unit for 3D arrays
-SPS5D = 1/20; % Sampling per time unit for 5D arrays
+SPS3D = 1/2; % Sampling per time unit for 3D arrays
+SPS5D = 1/200; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints/10
JOB2LOAD= -1;
%% OPTIONS AND NAMING
@@ -30,9 +37,9 @@ JOB2LOAD= -1;
% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; 4 : Coulomb; +/- for GK/DK)
CO = 1;
CLOS = 0; % Closure model (0: =0 truncation)
-NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
-% SIMID = 'nonlin_FCGK'; % Name of the simulation
-SIMID = 'test_even_p'; % Name of the simulation
+NL_CLOS = 0; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
+% SIMID = 'Linear_Device'; % Name of the simulation
+SIMID = 'simulation_A'; % Name of the simulation
% SIMID = ['v3.0_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation
NON_LIN = 1; % activate non-linearity (is cancelled if KXEQ0 = 1)
% INIT options
@@ -53,17 +60,14 @@ PMAXI = P; % Highest ion Hermite polynomial degree
JMAXI = J; % Highest '' Laguerre ''
KERN = 0; % Kernel model (0 : GK)
KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
-KXEQ0 = 0; % put kx = 0
KPAR = 0.0; % Parellel wave vector component
LAMBDAD = 0.0;
-kmax = N*pi/L;% Highest fourier mode
+kmax = Nx*pi/L;% Highest fourier mode
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
NOISE0 = 1.0e-5;
BCKGD0 = 0.0; % Init background
TAU = 1.0; % e/i temperature ratio
-ETAT = 0.0; % Temperature gradient
-ETAB = 1.0; % Magnetic gradient (1.0 to set R=LB)
INIT_PHI= 1; % Start simulation with a noisy phi and moments
MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 64959f91fad0e6143610d831fa34d532c8d16e10..15674da917cff5b59bde7d7b7e950bc957218d11 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -1,49 +1,50 @@
clear all;
addpath(genpath('../matlab')) % ... add
SUBMIT = 1; % To submit the job automatically
-CHAIN = 1; % To chain jobs (CHAIN = n will launch n jobs in chain)
+CHAIN = 2; % To chain jobs (CHAIN = n will launch n jobs in chain)
% EXECNAME = 'helaz_dbg';
- EXECNAME = 'helaz_3.81';
-for ETAN = [1/0.6]
+ EXECNAME = 'helaz_3.9';
+for K_N = [1/0.6]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CLUSTER PARAMETERS
CLUSTER.PART = 'prod'; % dbg or prod
% CLUSTER.PART = 'dbg';
-CLUSTER.TIME = '20:00:00'; % allocation time hh:mm:ss
+CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss
if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME = '00:30:00'; end;
CLUSTER.MEM = '128GB'; % Memory
CLUSTER.JNAME = 'HeLaZ';% Job name
-NP_P = 1; % MPI processes along p
-NP_KX = 48; % MPI processes along kx
+NP_P = 2; % MPI processes along p
+NP_KX = 24; % MPI processes along kx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
-ETAN = 1.0/0.6; % Density gradient drive (R/Ln)
+K_N = 1.0/0.6; % Density gradient drive (R/Ln)
NU_HYP = 0.0;
%% GRID PARAMETERS
-N = 500; % Frequency gridpoints (Nkx = N/2)
-L = 120; % Size of the squared frequency domain
+Nx = 300; % Realspace x-gridpoints
+Ny = 300; % Realspace y-gridpoints
+Lx = 120; % Size of the squared frequency domain
+Ly = 120; % Size of the squared frequency domain
Nz = 1; % number of perpendicular planes (parallel grid)
q0 = 1.0; % q factor ()
shear = 0.0; % magnetic shear
eps = 0.0; % inverse aspect ratio
-P = 4;
-J = 2;
+P = 8;
+J = 4;
%% TIME PARAMETERS
TMAX = 10000; % Maximal time unit
-DT = 1e-3; % Time step
+DT = 8e-3; % Time step
SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
SPS3D = 1/2; % Sampling per time unit for 3D arrays
SPS5D = 1/100; % Sampling per time unit for 5D arrays
-SPSCP = 0; % Sampling per time unit for checkpoints/10
-JOB2LOAD= 0; % start from t=0 if <0, else restart from outputs_$job2load
+JOB2LOAD= 2; % start from t=0 if <0, else restart from outputs_$job2load
%% OPTIONS AND NAMING
% Collision operator
% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; +/- for GK/DK)
-CO = 1;
+CO = 3;
CLOS = 0; % Closure model (0: =0 truncation)
NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
% SIMID = 'test_chained_job'; % Name of the simulation
@@ -56,7 +57,7 @@ INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0;
%% OUTPUTS
W_DOUBLE = 1;
W_GAMMA = 1; W_HF = 1;
-W_PHI = 1; W_NA00 = 1;
+W_PHI = 1; W_NA00 = 1;
W_DENS = 1; W_TEMP = 1;
W_NAPJ = 1; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -71,15 +72,14 @@ KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
KXEQ0 = 0; % put kx = 0
KPAR = 0.0; % Parellel wave vector component
LAMBDAD = 0.0;
-kmax = N*pi/L;% Highest fourier mode
+kmax = Nx*pi/Lx;% Highest fourier mode
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
% kmaxcut = 2.5;
MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
NOISE0 = 1.0e-5;
BCKGD0 = 0.0; % Init background
TAU = 1.0; % e/i temperature ratio
-ETAT = 0.0; % Temperature gradient
-ETAB = 1.0; % Magnetic gradient (1.0 to set R=LB)
+K_T = 0.0; % Temperature gradient
INIT_PHI= 1; % Start simulation with a noisy phi and moments
MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
diff --git a/wk/plot_cosol_mat.m b/wk/plot_cosol_mat.m
index bade9a4d4e2059b2db6671e0941f769b323eda75..6c5dcf6775c7fcfd89d26cb6a504748dbe4ce6f1 100644
--- a/wk/plot_cosol_mat.m
+++ b/wk/plot_cosol_mat.m
@@ -92,12 +92,9 @@ subplot(224)
%% Single eigenvalue analysis
-% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_6_3_NFLR_8_kp_4_bjf/scanfiles_00000/self.0.h5';
-% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_4_2_NFLR_2_kp_4/scanfiles_00000/self.0.h5';
-mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_6_P_4_J_2_N_75_kpm_6.0/scanfiles_00049/self.0.h5';
-% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0/scanfiles_00049/self.0.h5';
-% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_16_P_4_J_2_N_50_kpm_4.0/scanfiles_00049/self.0.h5';
-% mat_file_name = '/home/ahoffman/cosolver/trunc_gk.coulomb_NFLR_4_P_4_J_2_N_50_kpm_4.0/scanfiles_00049/self.0.h5';
+% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_20_P_4_J_2_N_50_kpm_4.0/scanfiles_00005/self.0.h5';
+mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_20_P_6_J_3_N_50_kpm_4.0/scanfiles_00036self.0.h5';
+
matidx = 01;
matidx = sprintf('%5.5i',matidx);disp(matidx);
@@ -121,7 +118,6 @@ mfns = {'/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
'/home/ahoffman/cosolver/gk.coulomb_NFLR_6_P_4_J_2_N_75_kpm_6.0.h5',...
'/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5',...
- '/home/ahoffman/cosolver/gk.coulomb_NFLR_16_P_4_J_2_N_50_kpm_4.0.h5',...
};
CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'FC 4 2 NFLR 16'};
figure
diff --git a/wk/plot_phi_ballooning.m b/wk/plot_phi_ballooning.m
new file mode 100644
index 0000000000000000000000000000000000000000..24011c8285b3b456296a9fc4d9da112ccf2ad8d5
--- /dev/null
+++ b/wk/plot_phi_ballooning.m
@@ -0,0 +1,65 @@
+phi_real=(real(PHI(:,:,:,end)));
+phi_imag=(imag(PHI(:,:,:,end)));
+% Apply baollooning tranform
+for iky=2
+ dims = size(phi_real);
+ phib_real = zeros( dims(1)*Nz ,1);
+ phib_imag= phib_real;
+ b_angle = phib_real;
+
+ midpoint = floor((dims(1)*Nz )/2)+1;
+ for ip =1: dims(1)
+ start_ = (ip -1)*Nz +1;
+ end_ = ip*Nz;
+ phib_real(start_:end_) = phi_real(ip,iky,:);
+ phib_imag(start_:end_) = phi_imag(ip,iky, :);
+ end
+
+ % Define ballooning angle
+ Nkx = numel(kx)-1; coordz = z;
+ idx = -Nkx:1:Nkx;
+ for ip =1: dims(1)
+ for iz=1:Nz
+ ii = Nz*(ip -1) + iz;
+ b_angle(ii) =coordz(iz) + 2*pi*idx(ip);
+ end
+ end
+
+ % normalize real and imaginary parts at chi =0
+ [~,idxLFS] = min(abs(b_angle -0));
+ phib = phib_real(:) + 1i * phib_imag(:);
+ % Normalize to the outermid plane
+ phib_norm = phib / phib( idxLFS) ;
+ phib_real_norm = real( phib_norm);%phib_real(:)/phib_real(idxLFS);
+ phib_imag_norm = imag( phib_norm);%phib_imag(:)/ phib_imag(idxLFS);
+
+
+ figure; hold all;
+ plot(b_angle/pi, phib_real_norm,'-b');
+ plot(b_angle/pi, phib_imag_norm ,'-r');
+ plot(b_angle/pi, phib_real_norm .^2 + phib_imag_norm.^2,'-k');
+ xlabel('$\chi / \pi$')
+ ylabel('$\phi_B (\chi)$');
+ title(['$(P,J) =(',num2str(PMAXI),', ', num2str(JMAXI),'$), $\nu =',num2str(NU),...
+ '$, $\epsilon = ',num2str(eps),'$, $k_y = ', num2str(ky( iky)),'$, $q =',num2str(q0),'$, $s = ', num2str(shear),'$, $R_N = ', ...
+ num2str(K_N),'$, $R_{T_i} = ', num2str(K_T),'$, $N_z =',num2str(Nz),'$']);
+ %set(gca,'Yscale','log')
+ %
+
+% %Check symmetry of the mode at the outter mid plane
+% figure; hold all;
+% right = phib_real(midpoint+1:end);
+% left = fliplr(phib_real(1:midpoint-1)');
+% up_down_symm = right(1:end) - left(1:end-1)';
+% %plot(b_angle(midpoint+1:end)/pi,phib_real(midpoint+1:end),'-xb');
+% plot(b_angle(midpoint+1:end)/pi,up_down_symm ,'-xb');
+ %plot(abs(b_angle(1:midpoint-1))/pi,phib_real(1:midpoint-1),'-xb');
+ %
+ %
+ % figure; hold all
+ % plot(b_angle/pi, phib_imag.^2 + phib_real.^2 ,'xk');
+ % %set(gca,'Yscale','log')
+ % xlabel('$\chi / \pi$')
+ % ylabel('$\phi_B (\chi)$');
+ % title(['$(P,J) =(',num2str(pmax),', ', num2str(jmax),'$), $\nu =',num2str(nu),'$, $\epsilon = ',num2str(epsilon),'$, $q =',num2str(safety_fac),'$, $s = ', num2str(shear),'$, $k_y =',num2str(ky),'$']);
+end
diff --git a/wk/shearless_linear_fluxtube.m b/wk/shearless_linear_fluxtube.m
new file mode 100644
index 0000000000000000000000000000000000000000..70311b3b5544f89ce5a41d8ab1c7d99862b1c144
--- /dev/null
+++ b/wk/shearless_linear_fluxtube.m
@@ -0,0 +1,113 @@
+RUN = 1; % To run or just to load
+addpath(genpath('../matlab')) % ... add
+default_plots_options
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set Up parameters
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% PHYSICAL PARAMETERS
+NU = 0.0; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_N = 2.22; % Density gradient drive
+K_T = 6.0; % Temperature '''
+SIGMA_E = 0.05196; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+%% GRID PARAMETERS
+Nx = 1; % real space x-gridpoints
+Ny = 2; % '' y-gridpoints
+Lx = 0; % Size of the squared frequency domain
+Ly = 2*pi/0.25; % Size of the squared frequency domain
+Nz = 24; % number of perpendicular planes (parallel grid)
+q0 = 2.7; % safety factor
+shear = 0.0; % magnetic shear
+eps = 0.18; % inverse aspect ratio
+%% TIME PARMETERS
+TMAX = 5; % Maximal time unit
+DT = 1e-3; % Time step
+SPS0D = 1; % Sampling per time unit for 2D arrays
+SPS2D = 0; % Sampling per time unit for 2D arrays
+SPS3D = 1; % Sampling per time unit for 2D arrays
+SPS5D = 1/100; % Sampling per time unit for 5D arrays
+SPSCP = 0; % Sampling per time unit for checkpoints
+JOB2LOAD= -1;
+%% OPTIONS
+SIMID = 'shearless_fluxtube'; % Name of the simulation
+% Collision operator
+% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle, 4 : Full Couloumb ; +/- for GK/DK)
+CO = 1;
+INIT_ZF = 0; ZF_AMP = 0.0;
+CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))
+NL_CLOS =-1; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
+KERN = 0; % Kernel model (0 : GK)
+%% OUTPUTS
+W_DOUBLE = 0;
+W_GAMMA = 1; W_HF = 1;
+W_PHI = 1; W_NA00 = 1;
+W_DENS = 1; W_TEMP = 1;
+W_NAPJ = 1; W_SAPJ = 0;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% unused
+HD_CO = 0.5; % Hyper diffusivity cutoff ratio
+kmax = Nx*pi/Lx;% Highest fourier mode
+NU_HYP = 0.0; % Hyperdiffusivity coefficient
+MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
+MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
+MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
+K_E = 0.00; % Electrostat '''
+GRADB = 1.0;
+CURVB = 1.0;
+INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0; INIT_PHI= 0;
+NOISE0 = 0.0e-4; % Init noise amplitude
+BCKGD0 = 1.0e-4; % Init background
+LAMBDAD = 0.0;
+KXEQ0 = 0; % put kx = 0
+NON_LIN = 0; % activate non-linearity (is cancelled if KXEQ0 = 1)
+%% PARAMETER SCANS
+
+if 1
+%% Parameter scan over PJ
+% PA = [2 4];
+% JA = [1 2];
+PA = [4];
+JA = [2];
+DTA= DT*ones(size(JA));%./sqrt(JA);
+% DTA= DT;
+mup_ = MU_P;
+muj_ = MU_J;
+Nparam = numel(PA);
+param_name = 'PJ';
+gamma_Ni00 = zeros(Nparam,floor(Nx/2)+1);
+gamma_Nipj = zeros(Nparam,floor(Nx/2)+1);
+gamma_phi = zeros(Nparam,floor(Nx/2)+1);
+% Ni00_ST = zeros(Nparam,floor(Nx/2)+1,TMAX/SPS3D);
+% PHI_ST = zeros(Nparam,floor(Nx/2)+1,TMAX/SPS3D);
+for i = 1:Nparam
+ % Change scan parameter
+ PMAXE = PA(i); PMAXI = PA(i);
+ JMAXE = JA(i); JMAXI = JA(i);
+ DT = DTA(i);
+ setup
+ system(['rm fort*.90']);
+ % Run linear simulation
+ if RUN
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk'])
+ end
+% Load and process results
+ %%
+ filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5'];
+ load_results
+end
+
+if 1
+%% Plot
+plot_phi_ballooning
+end
+end
+
+if 0
+%% Trajectories of some modes
+figure;
+for i = 1:10:Nx/2+1
+ semilogy(Ts3D,squeeze(abs(Ne00(i,2,1,:))),'DisplayName',['k=',num2str(kx(i))]); hold on;
+end
+end
\ No newline at end of file