diff --git a/matlab/LinearFit_s.m b/matlab/LinearFit_s.m
index 0d683f4f027e5a34c9232e92ecf21508153b467f..dc424911321e32633b28ae9fe48907c43ab3946b 100644
--- a/matlab/LinearFit_s.m
+++ b/matlab/LinearFit_s.m
@@ -1,4 +1,4 @@
-function [gamma,fit] = LinearFit_s(time,Na00abs, tstart, tend)
+function [gamma,fit] = LinearFit_s(time,Na00abs)
% LinearFit computes the growth rate and frequency from the time evolution of Napj
% - adapted from MOLI (B.J. Frei)
%
@@ -6,34 +6,29 @@ function [gamma,fit] = LinearFit_s(time,Na00abs, tstart, tend)
% ... amplitude ratio method
% We compute the mean of the growth rate over a time window
- Trun = time(end);
- if tstart == -1
- [~,Tstart_ind] = min(abs(time - 0.5*Trun));
- else
- [~,Tstart_ind] = min(abs(time - max(tstart,time(1))));
- end
-
- if tend == -1
- Tend_ind = numel(time)-1;
- else
- [~,Tend_ind] = min(abs(time - min(tend,Trun)));
- end
-
+% Trun = time(end);
+% if tstart == -1
+% [~,Tstart_ind] = min(abs(time - 0.5*Trun));
+% else
+% [~,Tstart_ind] = min(abs(time - max(tstart,time(1))));
+% end
+%
+% if tend == -1
+% Tend_ind = numel(time)-1;
+% else
+% [~,Tend_ind] = min(abs(time - min(tend,Trun)));
+% end
+%
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
+% gamma = mean(gammaoft(Tstart_ind:Tend_ind-1)); % ... take the mean of gamma over the time window
+ gamma = mean(gammaoft); % ... take the mean of gamma over the time window
% Return gamma(t) for amplitude ratio method
fit.gammaoft = gammaoft;
-
- % Return fit
- fit.t_min = tstart;
- fit.t_max = tend;
- fit.it_min = Tstart_ind;
- fit.it_max = Tend_ind;
end % ... end function
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index a2689bc1b83a42fcc22106895dd80d3ef71a8fc5..a8d4f665409ffe2911b5ae9fa84446bb5d83552f 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -2,9 +2,10 @@ CONTINUE = 1;
JOBNUM = 0; JOBFOUND = 0;
TJOB_SE = []; % Start and end times of jobs
NU_EVOL = []; % evolution of parameter nu between jobs
+MU_EVOL = []; % evolution of parameter nu between jobs
ETAB_EVOL= []; %
L_EVOL = []; %
-
+DT_EVOL = []; %
% FIELDS
Nipj_ = []; Nepj_ = [];
Ni00_ = []; Ne00_ = [];
@@ -18,8 +19,8 @@ Sipj_ = []; Sepj_ = [];
Pe_old = 1e9; Pi_old = Pe_old; Je_old = Pe_old; Ji_old = Pe_old;
while(CONTINUE)
- filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM);
- if exist(filename, 'file') == 2
+ filename = sprintf([BASIC.MISCDIR,'outputs_%.2d.h5'],JOBNUM);
+ if (exist(filename, 'file') == 2 && JOBNUM <= JOBNUMMAX)
% Load results of simulation #JOBNUM
load_results
% Check polynomials degrees
@@ -78,19 +79,20 @@ while(CONTINUE)
load_params
TJOB_SE = [TJOB_SE Ts0D(1) Ts0D(end)];
NU_EVOL = [NU_EVOL NU NU];
+ MU_EVOL = [MU_EVOL MU MU];
ETAB_EVOL = [ETAB_EVOL ETAB ETAB];
L_EVOL = [L_EVOL L L];
+ DT_EVOL = [DT_EVOL DT_SIM DT_SIM];
JOBFOUND = JOBFOUND + 1;
LASTJOB = JOBNUM;
- elseif (JOBNUM > 20)
+ Pe_old = Pe_new; Je_old = Je_new;
+ Pi_old = Pi_new; Ji_old = Ji_new;
+ elseif (JOBNUM > JOBNUMMAX)
CONTINUE = 0;
disp(['found ',num2str(JOBFOUND),' results']);
end
JOBNUM = JOBNUM + 1;
- Pe_old = Pe_new; Je_old = Je_new;
- Pi_old = Pi_new; Ji_old = Ji_new;
-
end
GGAMMA_RI = GGAMMA_; PGAMMA_RI = PGAMMA_; Ts0D = Ts0D_;
Nipj = Nipj_; Nepj = Nepj_; Ts5D = Ts5D_;
diff --git a/matlab/create_gif_5D.m b/matlab/create_gif_5D.m
index d77fd1f6c61c2ee467fffc56e97b68e0dbfad717..3acb6293da6814815df238ec07e268d0ac7c21e8 100644
--- a/matlab/create_gif_5D.m
+++ b/matlab/create_gif_5D.m
@@ -16,7 +16,7 @@ fig = figure('Color','white','Position', [100, 100, sz(2)*400, 400]);
else
yticks([])
end
- LEGEND = ['ln$|',FIELDNAME,'^{p',num2str(ij_-1),'}|$']; title(LEGEND);
+ LEGEND = ['$|',FIELDNAME,'^{p',num2str(ij_-1),'}|$']; title(LEGEND);
end
% colormap gray
axis tight manual % this ensures that getframe() returns a consistent size
@@ -35,7 +35,7 @@ fig = figure('Color','white','Position', [100, 100, sz(2)*400, 400]);
else
yticks([])
end
- LEGEND = ['ln$|',FIELDNAME,'^{p',num2str(ij_-1),'}|$'];
+ LEGEND = ['$|',FIELDNAME,'^{p',num2str(ij_-1),'}|$'];
title([LEGEND,', amp = ',sprintf('%.1e',scale)]);
end
suptitle(['$t \approx$', sprintf('%.3d',ceil(T(n)))]);
diff --git a/matlab/create_mov.m b/matlab/create_mov.m
new file mode 100644
index 0000000000000000000000000000000000000000..66da1aac2742e96ea68587c0cb89529a526b9d94
--- /dev/null
+++ b/matlab/create_mov.m
@@ -0,0 +1,51 @@
+title1 = GIFNAME;
+FIGDIR = BASIC.RESDIR;
+GIFNAME = [FIGDIR, GIFNAME,'.mp4'];
+vidfile = VideoWriter(GIFNAME,'Uncompressed AVI');
+open(vidfile);
+% Set colormap boundaries
+hmax = max(max(max(FIELD)));
+hmin = min(min(min(FIELD)));
+
+flag = 0;
+if hmax == hmin
+ disp('Warning : h = hmin = hmax = const')
+else
+% Setup figure frame
+figure('Color','white','Position', [100, 100, 400, 400]);
+ pcolor(X,Y,FIELD(:,:,1)); % to set up
+ colormap gray
+ axis tight manual % this ensures that getframe() returns a consistent size
+ if INTERP
+ shading interp;
+ end
+ in = 1;
+ nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:)));
+ for n = FRAMES % loop over selected frames
+ scale = max(max(abs(FIELD(:,:,n))));
+ pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot
+ if INTERP
+ shading interp;
+ 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)]);
+ drawnow
+ % Capture the plot as an image
+ frame = getframe(gcf);
+ writeVideo(vidfile,frame);
+ % terminal info
+ while nbytes > 0
+ fprintf('\b')
+ nbytes = nbytes - 1;
+ end
+ nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:)));
+ in = in + 1;
+ end
+ disp(' ')
+ disp(['Video saved @ : ',GIFNAME])
+ close(vidfile);
+end
+
diff --git a/matlab/load_marconi.m b/matlab/load_marconi.m
index 9540a049a13d66f2181e1a29ac44a2704414232e..fa714cbb18ae2f2ac69654c847f1d3dc9143bcf3 100644
--- a/matlab/load_marconi.m
+++ b/matlab/load_marconi.m
@@ -3,9 +3,14 @@ function [ RESDIR ] = load_marconi( outfilename )
% Detailed explanation goes here
hostfolder = outfilename(1:end-8);
hostfile = [hostfolder,'/out*'];
+ MISCDIR = ['/misc/HeLaZ_outputs/results/',outfilename(46:end-8),'/'];
RESDIR = ['../',outfilename(46:end-8),'/'];
+ miscfolder = [MISCDIR,'.'];
localfolder = [RESDIR,'.'];
- CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfile,' ',localfolder];
+ CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfile,' ',miscfolder];
+ disp(CMD);
+ system(CMD);
+ CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfolder,'/fort.90',' ',miscfolder];
disp(CMD);
system(CMD);
end
diff --git a/matlab/load_params.m b/matlab/load_params.m
index 7199be5b29a04543d0c53f6148cba3d8eaf70483..5443658e8856e51f84100e3ce2e0929fa49a34a6 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -12,7 +12,8 @@ NR = h5readatt(filename,'/data/input','nr');
NZ = h5readatt(filename,'/data/input','nz');
L = h5readatt(filename,'/data/input','Lr');
CLOS = h5readatt(filename,'/data/input','CLOS');
-MU = str2num(filename(end-18:end-14));
+DT_SIM = h5readatt(filename,'/data/input','dt');
+MU = str2num(filename(end-18:end-14));
W_GAMMA = h5readatt(filename,'/data/input','write_gamma') == 'y';
W_PHI = h5readatt(filename,'/data/input','write_phi') == 'y';
W_NA00 = h5readatt(filename,'/data/input','write_Na00') == 'y';
diff --git a/matlab/load_results.m b/matlab/load_results.m
index fc41f9c4f13a013dc996d153051a3ca1d61c0dac..c4489bc0c1c7f62f0ac015c575b5a841b287a594 100644
--- a/matlab/load_results.m
+++ b/matlab/load_results.m
@@ -1,5 +1,4 @@
%% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM);
disp(['Loading ',filename])
% Loading from output file
CPUTIME = h5readatt(filename,'/data/input','cpu_time');
@@ -7,6 +6,13 @@ DT_SIM = h5readatt(filename,'/data/input','dt');
[Pe, Je, Pi, Ji, kr, kz] = load_grid_data(filename);
+W_GAMMA = strcmp(h5readatt(filename,'/data/input','write_gamma'),'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');
+W_SAPJ = strcmp(h5readatt(filename,'/data/input','write_Sapj') ,'y');
+
+
if W_GAMMA
[ GGAMMA_RI, Ts0D, dt0D] = load_0D_data(filename, 'gflux_ri');
PGAMMA_RI = load_0D_data(filename, 'pflux_ri');
diff --git a/matlab/photomaton.m b/matlab/photomaton.m
new file mode 100644
index 0000000000000000000000000000000000000000..a73cc5f8c244879ec8d7c78631bcb9a9a3c33d17
--- /dev/null
+++ b/matlab/photomaton.m
@@ -0,0 +1,67 @@
+if 0
+%% Photomaton : real space
+% FIELD = ni00; FNAME = 'ni'; XX = RR; YY = ZZ;
+% FIELD = phi; FNAME = 'phi'; XX = RR; YY = ZZ;
+FIELD = dr2phi; FNAME = 'dr2phi'; XX = RR; YY = ZZ;
+% FIELD = fftshift(abs(Ni00),2); FNAME = 'Fni'; XX = fftshift(KR,2); YY = fftshift(KZ,2);
+% FIELD = fftshift(abs(PHI),2); FNAME = 'Fphi'; XX = fftshift(KR,2); YY = fftshift(KZ,2);
+tf = 150; [~,it1] = min(abs(Ts2D-tf));
+tf = 300; [~,it2] = min(abs(Ts2D-tf));
+tf = 800; [~,it3] = min(abs(Ts2D-tf));
+tf = 1200; [~,it4] = min(abs(Ts2D-tf));
+fig = figure; FIGNAME = [FNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 350])
+plt = @(x) x;%./max(max(x));
+ subplot(141)
+ DATA = plt(FIELD(:,:,it1));
+ pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ colormap gray
+ xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$');set(gca,'ytick',[]);
+ title(sprintf('$t c_s/R=%.0f$',Ts2D(it1)));
+ subplot(142)
+ DATA = plt(FIELD(:,:,it2));
+ pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ colormap gray
+ xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]);
+ title(sprintf('$t c_s/R=%.0f$',Ts2D(it2)));
+ subplot(143)
+ DATA = plt(FIELD(:,:,it3));
+ pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ colormap gray
+ xlabel('$r/\rho_s$');ylabel('$z/\rho_s$');set(gca,'ytick',[]);
+ title(sprintf('$t c_s/R=%.0f$',Ts2D(it3)));
+ subplot(144)
+ DATA = plt(FIELD(:,:,it4));
+ pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ colormap gray
+ xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]);
+ title(sprintf('$t c_s/R=%.0f$',Ts2D(it4)));
+suptitle(['$\',FNAME,'$, $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
+ ', $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']);
+save_figure
+end
+
+%%
+if 0
+%% Show frame in kspace
+tf = 1000; [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf));
+fig = figure; FIGNAME = ['krkz_',sprintf('t=%.0f',Ts2D(it2)),'_',PARAMS];set(gcf, 'Position', [100, 100, 700, 600])
+CLIM = [0,1e3]
+ subplot(223); plt = @(x) fftshift((abs(x)),2);
+ pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(PHI(:,:,it2))); set(pclr, 'edgecolor','none');
+ caxis(CLIM); colormap hot
+ xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('$t c_s/R=%.0f$',Ts2D(it2))); legend('$|\hat\phi|$');
+ subplot(222); plt = @(x) fftshift(abs(x),2);
+ pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ni00(:,:,it2))); set(pclr, 'edgecolor','none');
+ caxis(CLIM); colormap hot
+ xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_i^{00}|$');
+ subplot(221); plt = @(x) fftshift(abs(x),2);
+ pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ne00(:,:,it2))); set(pclr, 'edgecolor','none');
+ caxis(CLIM); colormap hot
+ xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_e^{00}|$');
+ subplot(224); plt = @(x) fftshift((abs(x)),2);
+ colorbar;
+ caxis(CLIM); colormap hot
+% pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Si00(:,:,it5))); set(pclr, 'edgecolor','none'); colorbar;
+% xlabel('$k_r$'); ylabel('$k_z$');legend('$\hat S_i^{00}$');
+save_figure
+end
\ No newline at end of file
diff --git a/matlab/write_sbash_marconi.m b/matlab/write_sbash_marconi.m
index b4fdf01b51f22ea571a9eea7a20a38e200c3155c..e544d84c83af6f67b157847de6d3e3dcc15bffdf 100644
--- a/matlab/write_sbash_marconi.m
+++ b/matlab/write_sbash_marconi.m
@@ -14,7 +14,7 @@ fprintf(fid,[...
'cp $HOME/HeLaZ/wk/batch_script.sh .\n',...
...
'sbatch batch_script.sh\n',...
-'echo tail -f $CINECA_SCRATCH/HeLaZ/results/',BASIC.SIMID,'/',BASIC.PARAMS,'/out.txt']);
+'echo tail -f $CINECA_SCRATCH/HeLaZ',BASIC.RESDIR(3:end),'out.txt']);
fclose(fid);
system(['cp setup_and_run.sh ',BASIC.RESDIR,'/.']);
@@ -41,4 +41,5 @@ fprintf(fid,[...
fclose(fid);
system(['cp batch_script.sh ',BASIC.RESDIR,'/.']);
-system('scp {fort.90,setup_and_run.sh,batch_script.sh} ahoffman@login.marconi.cineca.it:/marconi/home/userexternal/ahoffman/HeLaZ/wk');
\ No newline at end of file
+system('scp {fort.90,setup_and_run.sh,batch_script.sh} ahoffman@login.marconi.cineca.it:/marconi/home/userexternal/ahoffman/HeLaZ/wk > trash.txt');
+system('rm trash.txt');
\ No newline at end of file
diff --git a/wk/ZF_fourier_analysis.m b/wk/ZF_fourier_analysis.m
new file mode 100644
index 0000000000000000000000000000000000000000..c63ea03c1e6fa89cc6f72c0df99d7fa4b5b9d03b
--- /dev/null
+++ b/wk/ZF_fourier_analysis.m
@@ -0,0 +1,73 @@
+%% Zonal flow spectral analysis
+fig = figure; FIGNAME = ['zonal_flow_spectral_analysis_',PARAMS];
+tend = Ts0D(end); tstart = tend - TAVG;
+[~,its0D] = min(abs(Ts0D-tstart));
+[~,ite0D] = min(abs(Ts0D-tend));
+set(gcf, 'Position', [100, 100, 800, 400])
+ % Time series analysis (burst period and time frequencies spectrum)
+ subplot(121)
+ samplerate = Ts0D(2)-Ts0D(1);
+ Y = log(PGAMMA_RI(its0D:ite0D)*(2*pi/Nr/Nz)^2);
+ [n,~] = size(Y);
+ Yy= fft(y);
+ Pot = Yy .* conj(Yy) / n;
+ Freq = (samplerate / n * (1:n))';
+ Pot(1) = 0;
+ nmax = min(200,round(n/2));
+ [amax, itmax] = max(Pot);
+ plot((0:nmax-1) , Pot(1:nmax)/amax,'DisplayName','$\Gamma_r(\omega)$');hold on;
+ plot([itmax-1,itmax-1],[0,1],'--k', 'DisplayName',['$T_{per}\approx',num2str(round(2*pi/Freq(itmax))),'\rho_s/c_s$']);
+ legend('show'); grid on; box on; xlabel('Period number'); yticks([]);
+ title('ZF temporal spectrum')
+ % Space analysis (spatial period of ZF)
+ subplot(122)
+ Yy = mean(mean(1i*KR.*PHI(:,:,its2D:ite2D),3),2); hold on;
+ [n,~] = size(Yy);
+ Pot = Yy .* conj(Yy) / n;
+ [amax, ikmax] = max(Pot);
+ nmax = 20;
+ plot(0:nmax,Pot(1:nmax+1)/amax,'DisplayName','$\tilde V_E(k_r)$');
+ plot([ikmax-1,ikmax-1],[0,1],'--k', 'DisplayName',['$L_z=',num2str(2*pi/kr(ikmax)),'\rho_s$']);
+ grid on; box on;
+ title('ZF spatial spectrum')
+ xlabel('mode number'); yticks([]); legend('show')
+save_figure
+
+%% Shear and phi amp phase space
+fig = figure; FIGNAME = ['phi_shear_phase_space_',PARAMS];
+set(gcf, 'Position', [100, 100, 700, 500])
+t1 = Ts2D(end); t0 = t1 - TAVG;
+[~,its2D] = min(abs(Ts2D-t0)); [~,ite2D] = min(abs(Ts2D-t1));
+scatter(phi_maxr_avgz(its2D:ite2D),shear_maxr_avgz(its2D:ite2D),35,Ts2D(its2D:ite2D),'.',...
+ 'DisplayName',PARAMS); cbar = colorbar;ylabel(cbar,'$t c_s/\rho_s$','Interpreter','LaTeX')
+hold on
+xlabel('$\langle \phi \rangle_z^r$'); ylabel('$\langle dV_E/dr \rangle_z^r$')
+grid on; title('ES pot. vs Shear phase space')
+% plot(phi_avgr_maxz(its2D:ite2D),shear_avgr_maxz(its2D:ite2D),'-')
+% plot(phi_maxr_maxz(its2D:ite2D),shear_maxr_maxz(its2D:ite2D),'-')
+% plot(phi_avgr_avgz(its2D:ite2D),shear_avgr_avgz(its2D:ite2D),'-')
+save_figure
+
+%% Non zonal quantities
+PHI_NZ = PHI;
+PHI_NZ(ikmax-1:ikmax+1,:,:) = 0;
+
+phi_nz = zeros(Nr,Nz,Ns2D);
+for it = 1:numel(Ts2D)
+ PH_ = PHI_NZ(:,:,it);
+ phi_nz (:,:,it) = real(fftshift(ifft2((PH_),Nr,Nz)));
+end
+%%
+t0 = 1000;
+[~, it02D] = min(abs(Ts2D-t0));
+[~, it05D] = min(abs(Ts5D-t0));
+skip_ = 10;
+DELAY = 0.005*skip_;
+FRAMES_2D = it02D:skip_:numel(Ts2D);
+if 0
+%% Phi non zonal real space
+GIFNAME = ['phi_nz',sprintf('_%.2d',JOBNUM),'_',PARAMS];INTERP = 0;
+FIELD = real(phi_nz); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
+FIELDNAME = '$\phi_{NZ}$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
+create_gif
+end
\ No newline at end of file
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 4eca5eacb82d3c3e3ca973bbe852df520bdbd859..e66ebd2a12931ef09ed498f381e7c503e78c9eb3 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -1,23 +1,28 @@
+addpath(genpath('../matlab')) % ... add
%% Load results
outfile ='';
-if 0
- %% Load from Marconi
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.4_eta_0.6_nu_1e-01/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_80_P_10_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_4e-03/out.txt';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e+00/200x100_L_120_P_12_J_6_eta_0.6_nu_1e+00_DGGK_CLOS_0_mu_2e-02/out.txt';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_12_J_6_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_12_J_6_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_1e-02/out.txt';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_12_J_6_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
- BASIC.RESDIR = load_marconi(outfile);
+if 0 % Local results
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='';
+ BASIC.RESDIR = ['../results/',outfile,'/'];
+ BASIC.MISCDIR = ['../results/',outfile,'/'];
end
-if 0
- %% Load from Daint
- outfile ='/scratch/snx3000/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_12_J_6_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
- BASIC.RESDIR = load_daint(outfile);
+if 1 % Marconi results
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.9_nu_5e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
+% load_marconi(outfile);
+BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
+BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
end
+
%%
% JOBNUM = 1; load_results;
-compile_results
+JOBNUMMAX = 20; compile_results %Compile the results from first output found to JOBNUMMAX if existing
%% Retrieving max polynomial degree and sampling info
Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe);
@@ -91,9 +96,6 @@ end
% Post processing
disp('- post processing')
-E_pot = zeros(1,Ns2D); % Potential energy n^2
-E_kin = zeros(1,Ns2D); % Kinetic energy grad(phi)^2
-ExB = zeros(1,Ns2D); % ExB drift intensity \propto |\grad \phi|
% gyrocenter and particle flux from real space
GFlux_ri = zeros(1,Ns2D); % Gyrocenter flux Gamma = <ni drphi>
GFlux_zi = zeros(1,Ns2D); % Gyrocenter flux Gamma = <ni dzphi>
@@ -111,6 +113,12 @@ phi_maxr_maxz = zeros(1,Ns2D); % Time evol. of the norm of phi
phi_avgr_maxz = zeros(1,Ns2D); % Time evol. of the norm of phi
phi_maxr_avgz = zeros(1,Ns2D); % Time evol. of the norm of phi
phi_avgr_avgz = zeros(1,Ns2D); % Time evol. of the norm of phi
+
+shear_maxr_maxz = zeros(1,Ns2D); % Time evol. of the norm of shear
+shear_avgr_maxz = zeros(1,Ns2D); % Time evol. of the norm of shear
+shear_maxr_avgz = zeros(1,Ns2D); % Time evol. of the norm of shear
+shear_avgr_avgz = zeros(1,Ns2D); % Time evol. of the norm of shear
+
Ne_norm = zeros(Npe,Nje,Ns5D); % Time evol. of the norm of Napj
Ni_norm = zeros(Npi,Nji,Ns5D); % .
@@ -122,7 +130,12 @@ for it = 1:numel(Ts2D) % Loop over 2D arrays
phi_avgr_maxz(it) = max(mean(squeeze(phi(:,:,it))));
phi_maxr_avgz(it) = mean( max(squeeze(phi(:,:,it))));
phi_avgr_avgz(it) = mean(mean(squeeze(phi(:,:,it))));
- ExB(it) = max(max(max(abs(phi(3:end,:,it)-phi(1:end-2,:,it))/(2*dr))),max(max(abs(phi(:,3:end,it)-phi(:,1:end-2,it))'/(2*dz))));
+
+ shear_maxr_maxz(it) = max( max(squeeze(-(dr2phi(:,:,it)))));
+ shear_avgr_maxz(it) = max(mean(squeeze(-(dr2phi(:,:,it)))));
+ shear_maxr_avgz(it) = mean( max(squeeze(-(dr2phi(:,:,it)))));
+ shear_avgr_avgz(it) = mean(mean(squeeze(-(dr2phi(:,:,it)))));
+
GFlux_ri(it) = sum(sum(ni00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz;
GFlux_zi(it) = sum(sum(-ni00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz;
GFlux_re(it) = sum(sum(ne00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz;
@@ -145,15 +158,19 @@ disp('- growth rate')
[tmp,tmax] = max(GGAMMA_RI*(2*pi/Nr/Nz)^2);
[~,itmax] = min(abs(Ts2D-tmax));
tstart = 0.1 * Ts2D(itmax); tend = 0.5 * Ts2D(itmax);
+[~,its2D_lin] = min(abs(Ts2D-tstart));
+[~,ite2D_lin] = min(abs(Ts2D-tend));
+
g_ = zeros(Nkr,Nkz);
for ikr = 1:Nkr
for ikz = 1:Nkz
- g_(ikr,ikz) = LinearFit_s(Ts2D,squeeze(abs(Ni00(ikr,ikz,:))),tstart,tend);
+ [g_(ikr,ikz), ~] = LinearFit_s(Ts2D(its2D_lin:ite2D_lin),squeeze(abs(Ni00(ikr,ikz,its2D_lin:ite2D_lin))));
end
end
[gmax,ikzmax] = max(g_(1,:));
kzmax = abs(kz(ikzmax));
Bohm_transport = ETAB/ETAN*gmax/kzmax^2;
+
%% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
default_plots_options
disp('Plots')
@@ -202,7 +219,7 @@ subplot(111);
plot(kz,g_(1,:),'-','DisplayName','Linear growth rate'); hold on;
plot([max(kz)*2/3,max(kz)*2/3],[0,10],'--k', 'DisplayName','2/3 Orszag AA');
grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma R/c_s$'); legend('show');
-% ylim([0,max(g_(1,:))]); xlim([0,max(kz)]);
+ ylim([0,max(g_(1,:))]); xlim([0,max(kz)]);
subplot(224)
clr = line_colors(min(ip,numel(line_colors(:,1))),:);
lstyle = line_styles(min(ij,numel(line_styles)));
@@ -214,183 +231,110 @@ subplot(111);
save_figure
end
-if 0
-%% Particle fluxes
-SCALING = Nkr*dkr * Nkz*dkz;
-fig = figure; FIGNAME = ['gamma',sprintf('_%.2d',JOBNUM),'_',PARAMS];
-set(gcf, 'Position', [100, 100, 800, 300])
- semilogy(Ts2D,GFLUX_RI, 'color', line_colors(2,:)); hold on
- plot(Ts5D,PFLUX_RI,'.', 'color', line_colors(2,:)); hold on
- plot(Ts2D,SCALING*GFlux_ri, 'color', line_colors(1,:)); hold on
- plot(Ts5D,SCALING*PFlux_ri,'.', 'color', line_colors(1,:)); hold on
- xlabel('$tc_{s0}/\rho_s$'); ylabel('$\Gamma_r$'); grid on
- title(['$\eta=',num2str(ETAB),'\quad',...
- '\nu_{',CONAME,'}=',num2str(NU),'$', ' $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$'])
- legend('Gyro Flux','Particle flux', 'iFFT GFlux', 'iFFT PFlux')%'$\eta\gamma_{max}/k_{max}^2$')
-save_figure
-end
-
-if 0
+if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-fig = figure; FIGNAME = ['space_time_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
+TAVG = 1500; % Averaging time duration
+%Compute steady radial transport
+tend = Ts0D(end); tstart = tend - TAVG;
+[~,its0D] = min(abs(Ts0D-tstart));
+[~,ite0D] = min(abs(Ts0D-tend));
+SCALE = (2*pi/Nr/Nz)^2;
+gamma_infty_avg = mean(PGAMMA_RI(its0D:ite0D))*SCALE;
+gamma_infty_std = std (PGAMMA_RI(its0D:ite0D))*SCALE;
+% Compute steady shearing rate
+tend = Ts2D(end); tstart = tend - TAVG;
+[~,its2D] = min(abs(Ts2D-tstart));
+[~,ite2D] = min(abs(Ts2D-tend));
+shear_infty_avg = mean(shear_maxr_avgz(its2D:ite2D));
+shear_infty_std = std (shear_maxr_avgz(its2D:ite2D));
+% plots
+fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
subplot(311)
- semilogy(Ts2D,GFLUX_RI,'-'); hold on
- plot(Ts5D,PFLUX_RI,'.'); hold on
-% plot(Ts2D,Bohm_transport*ones(size(Ts2D)),'--'); hold on
- ylabel('$\Gamma_r$'); grid on
- title(['$\eta=',num2str(ETAB),'\quad',...
- '\nu_{',CONAME,'}=',num2str(NU),'$'])
- legend(['$P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$'],'Particle flux')%'$\eta\gamma_{max}/k_{max}^2$')
-% set(gca,'xticklabel',[])
+ plot(Ts0D,PGAMMA_RI*(2*pi/Nr/Nz)^2,'DisplayName','$\langle n_i d\phi/dz \rangle_z$'); 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('Transport')
+ ylim([0,gamma_infty_avg*5.0]); xlim([0,Ts0D(end)]);
+ title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
+ ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
+ ' $\mu_{hd}=$',num2str(MU)]);
+ legend('show');
+
subplot(312)
- yyaxis left
- plot(Ts2D,squeeze(max(max((phi)))))
- ylabel('$\max \phi$')
- yyaxis right
- plot(Ts2D,squeeze(mean(max(dr2phi))))
- ylabel('$s\sim\langle\partial_r^2\phi\rangle_z$'); grid on
-% set(gca,'xticklabel',[])
+ clr = line_colors(1,:);
+ lstyle = line_styles(1);
+% plot(Ts2D,phi_maxr_maxz,'DisplayName','$\max_{r,z}(\phi)$'); hold on;
+% plot(Ts2D,phi_maxr_avgz,'DisplayName','$\max_{r}\langle \phi\rangle_z$'); hold on;
+% plot(Ts2D,phi_avgr_maxz,'DisplayName','$\max_{z}\langle \phi\rangle_r$'); hold on;
+% plot(Ts2D,phi_avgr_avgz,'DisplayName','$\langle \phi\rangle_{r,z}$'); hold on;
+ plot(Ts2D,shear_maxr_maxz,'DisplayName','$\max_{r,z}(s_\phi)$'); hold on;
+ plot(Ts2D,shear_maxr_avgz,'DisplayName','$\max_{r}\langle s_\phi\rangle_z$'); hold on;
+ plot(Ts2D,shear_avgr_maxz,'DisplayName','$\max_{z}\langle s_\phi\rangle_r$'); hold on;
+ plot(Ts2D,shear_avgr_avgz,'DisplayName','$\langle s_\phi\rangle_{r,z}$'); hold on;
+ plot(Ts2D(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([0,Ts0D(end)]);
+ grid on; ylabel('Shear amp.'); legend('show');set(gca,'xticklabel',[])
subplot(313)
[TY,TX] = meshgrid(r,Ts2D);
- pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,:),2))'); set(pclr, 'edgecolor','none'); %colorbar;
- xlabel('$t c_s/R$'), ylabel('$r/\rho_s$')
- legend('$\langle\partial_r \phi\rangle_z$')
+% pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,:),2))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_r\phi\rangle_z$') %colorbar;
+ pclr = pcolor(TX,TY,squeeze(mean(dr2phi(:,:,:),2))'); set(pclr, 'edgecolor','none'); legend('Shear ($\langle \partial_r^2\phi\rangle_z$)') %colorbar;
+ caxis(1*shear_infty_avg*[-1 1]); xlabel('$t c_s/R$'), ylabel('$r/\rho_s$'); colormap hot
save_figure
end
if 0
-%% Photomaton : real space
-% FIELD = ni00; FNAME = 'ni'; XX = RR; YY = ZZ;
-FIELD = phi; FNAME = 'phi'; XX = RR; YY = ZZ;
-% FIELD = fftshift(abs(Ni00),2); FNAME = 'Fni'; XX = fftshift(KR,2); YY = fftshift(KZ,2);
-% FIELD = fftshift(abs(PHI),2); FNAME = 'Fphi'; XX = fftshift(KR,2); YY = fftshift(KZ,2);
-tf = 100; [~,it1] = min(abs(Ts2D-tf));
-tf = 118; [~,it2] = min(abs(Ts2D-tf));
-tf = 140; [~,it3] = min(abs(Ts2D-tf));
-tf = 300; [~,it4] = min(abs(Ts2D-tf));
-fig = figure; FIGNAME = [FNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 400])
-plt = @(x) x;%./max(max(x));
- subplot(141)
- DATA = plt(FIELD(:,:,it1));
- pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap gray
- xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$');set(gca,'ytick',[]);
- title(sprintf('$t c_s/R=%.0f$',Ts2D(it1)));
- subplot(142)
- DATA = plt(FIELD(:,:,it2));
- pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap gray
- xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]);
- title(sprintf('$t c_s/R=%.0f$',Ts2D(it2)));
- subplot(143)
- DATA = plt(FIELD(:,:,it3));
- pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap gray
- xlabel('$r/\rho_s$');ylabel('$z/\rho_s$');set(gca,'ytick',[]);
- title(sprintf('$t c_s/R=%.0f$',Ts2D(it3)));
- subplot(144)
- DATA = plt(FIELD(:,:,it4));
- pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap gray
- xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]);
- title(sprintf('$t c_s/R=%.0f$',Ts2D(it4)));
-% suptitle(['$\',FNAME,'$, $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
-% ', $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']);
-save_figure
-end
-
-%%
-if 0
-%% Show frame in kspace
-tf = 0; [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf));
-fig = figure; FIGNAME = ['krkz_',sprintf('t=%.0f',Ts2D(it2)),'_',PARAMS];set(gcf, 'Position', [100, 100, 700, 600])
- subplot(221); plt = @(x) fftshift((abs(x)),2);
- pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(PHI(:,:,it2))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('$t c_s/R=%.0f$',Ts2D(it2))); legend('$|\hat\phi|$');
- subplot(222); plt = @(x) fftshift(abs(x),2);
- pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ni00(:,:,it2))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_i^{00}|$');
- subplot(223); plt = @(x) fftshift((abs(x)),2); FIELD = squeeze(Nipj(1,2,:,:,:));
- pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(FIELD(:,:,it5))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_i^{pj=01}|$');
- subplot(224); plt = @(x) fftshift((abs(x)),2);
- pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Si00(:,:,it5))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$k_r$'); ylabel('$k_z$');legend('$\hat S_i^{00}$');
-save_figure
-end
-
-%%
-if 0
-%% Hermite energy spectra
-% tf = Ts2D(end-3);
-fig = figure; FIGNAME = ['hermite_spectrum_',PARAMS];set(gcf, 'Position', [100, 100, 1800, 600]);
-plt = @(x) squeeze(x);
-for ij = 1:Nji
- subplotnum = 100+Nji*10+ij;
- subplot(subplotnum)
- for it5 = 1:2:Ns5D
- alpha = it5*1.0/Ns5D;
- loglog(Pi(1:2:end),plt(epsilon_i_pj(1:2:end,ij,it5)),...
- 'color',(1-alpha)*[0.8500, 0.3250, 0.0980]+alpha*[0, 0.4470, 0.7410],...
- 'DisplayName',['t=',num2str(Ts5D(it5))]); hold on;
- end
- grid on;
- xlabel('$p$');
- TITLE = ['$\sum_{kr,kz} |N_i^{p',num2str(Ji(ij)),'}|^2$']; title(TITLE);
-end
-save_figure
-end
-%%
-if 0
-%% Laguerre energy spectra
-% tf = Ts2D(end-3);
-fig = figure; FIGNAME = ['laguerre_spectrum_',PARAMS];set(gcf, 'Position', [100, 100, 500, 400]);
-plt = @(x) squeeze(x);
-for it5 = 1:2:Ns5D
- alpha = it5*1.0/Ns5D;
- loglog(Ji,plt(max(epsilon_i_pj(:,:,it5),[],1)),...
- 'color',(1-alpha)*[0.8500, 0.3250, 0.0980]+alpha*[0, 0.4470, 0.7410],...
- 'DisplayName',['t=',num2str(Ts5D(it5))]); hold on;
-end
-grid on;
-xlabel('$j$');
-TITLE = ['$\max_p\sum_{kr,kz} |N_i^{pj}|^2$']; title(TITLE);
+%% Space time diagramms
+tstart = 0; tend = Ts2D(end);
+[~,itstart] = min(abs(Ts2D-tstart));
+[~,itend] = min(abs(Ts2D-tend));
+trange = itstart:itend;
+[TY,TX] = meshgrid(r,Ts2D(trange));
+fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
+ subplot(211)
+ pclr = pcolor(TX,TY,squeeze(mean(ni00(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
+ shading interp
+ colormap hot;
+ caxis([0.0,0.02]);
+ xticks([]); ylabel('$r/\rho_s$')
+ legend('$\Gamma_r(r)=\langle n_i^{GC}\partial_z\phi\rangle_z$')
+ title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
+ ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
+ ' $\mu_{hd}=$',num2str(MU)]);
+ subplot(212)
+ pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$t c_s/R$'), ylabel('$r/\rho_s$')
+ legend('$\langle \partial_r\phi\rangle_z$')
save_figure
end
%%
-no_AA = (2:floor(2*Nkr/3));
-tKHI = 100;
-[~,itKHI] = min(abs(Ts2D-tKHI));
-after_KHI = (itKHI:Ns2D);
-if 0
-%% Phi frequency space time diagram at kz=0
-fig = figure; FIGNAME = ['phi_freq_diag_',PARAMS];set(gcf, 'Position', [100, 100, 500, 400]);
- [TY,TX] = meshgrid(Ts2D(after_KHI),kr(no_AA));
- pclr = pcolor(TX,TY,(squeeze(abs(PHI(no_AA,1,(after_KHI)))))); set(pclr, 'edgecolor','none'); colorbar;
- ylabel('$t c_s/R$'), xlabel('$0<k_r<2/3 k_r^{\max}$')
- legend('$\log|\tilde\phi(k_z=0)|$')
- title('Spectrogram of $\phi$')
-end
-%%
-t0 = 0;
+t0 = 0000;
[~, it02D] = min(abs(Ts2D-t0));
[~, it05D] = min(abs(Ts5D-t0));
-skip_ = 2;
-DELAY = 0.005*skip_;
+skip_ = 10;
+DELAY = 1e-3*skip_;
FRAMES_2D = it02D:skip_:numel(Ts2D);
FRAMES_5D = it05D:skip_:numel(Ts5D);
%% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 0
%% Density ion
-GIFNAME = ['ni',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 1;
+GIFNAME = ['ni',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
FIELDNAME = '$n_i$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
create_gif
end
if 0
+%% Density electrons
+GIFNAME = ['ne',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
+FIELD = real(ne00); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
+FIELDNAME = '$n_e$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
+% create_gif
+create_mov
+end
+if 0
%% Phi real space
-GIFNAME = ['phi',sprintf('_%.2d',JOBNUM),'_',PARAMS];INTERP = 1;
+GIFNAME = ['phi',sprintf('_%.2d',JOBNUM),'_',PARAMS];INTERP = 0;
FIELD = real(phi); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
FIELDNAME = '$\phi$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
create_gif
@@ -422,7 +366,7 @@ if 0
GIFNAME = ['phi_kz0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING = 0;
FIELD =squeeze(log10(abs(PHI(no_AA,1,:)))); linestyle = '-.'; FRAMES = FRAMES_2D;
X = kr(no_AA); T = Ts2D; YMIN = -30; YMAX = 6; XMIN = min(kr); XMAX = max(kr);
-FIELDNAME = '$|\tilde\phi(k_z=0)|$'; XNAME = '$k_r\rho_s$';
+FIELDNAME = '$\log_{10}|\tilde\phi(k_z=0)|$'; XNAME = '$k_r\rho_s$';
create_gif_1D
end
if 0
@@ -444,23 +388,23 @@ if 0
GIFNAME = ['Sip0_kr',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
plt = @(x) squeeze(max((abs(x)),[],4));
FIELD =plt(Sipj(:,1,:,:,:)); X = kr'; Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
-FIELDNAME = '$N_i^{p0}$'; XNAME = '$k_{max}\rho_s$'; YNAME = '$P$';
+FIELDNAME = '$N_i^{p0}$'; XNAME = '$k_{max}\rho_s$'; YNAME = '$P$, $\sum_z$';
create_gif_imagesc
end
if 0
%% maxkz, kr vs p, for all Nipj over time
GIFNAME = ['Nipj_kr',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
-plt = @(x) squeeze(max((abs(x)),[],4));
+plt = @(x) squeeze(sum((abs(x)),4));
FIELD = plt(Nipj); X = kr'; Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
-FIELDNAME = 'N_i'; XNAME = '$k_r\rho_s$'; YNAME = '$P$, ${k_z}^{max}$';
+FIELDNAME = 'N_i'; XNAME = '$k_r\rho_s$'; YNAME = '$P$';
create_gif_5D
end
if 0
%% maxkr, kz vs p, for all Nipj over time
GIFNAME = ['Nipj_kz',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
-plt = @(x) fftshift(squeeze(max((abs(x)),[],3)),3);
+plt = @(x) fftshift(squeeze(sum((abs(x)),3)),3);
FIELD = plt(Nipj); X = sort(kz'); Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
-FIELDNAME = 'N_i'; XNAME = '$k_z\rho_s$'; YNAME = '$P$, ${k_r}^{max}$';
+FIELDNAME = 'N_i'; XNAME = '$k_z\rho_s$'; YNAME = '$P$, $\sum_r$';
create_gif_5D
end
%%
\ No newline at end of file
diff --git a/wk/analysis_stat_2D.m b/wk/analysis_stat_2D.m
index 07d4dcf9fc8b5e579ec0cacf6e88f8196f155ac1..cae39cf8a58e7212bdf2e5c0b615dd4407684a52 100644
--- a/wk/analysis_stat_2D.m
+++ b/wk/analysis_stat_2D.m
@@ -1,7 +1,7 @@
if 0
%% ni ne phase space for different time windows
fig = figure;
- t0 = 80; t1 = 100; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
+ t0 = 3500; t1 = 4500; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
plt = @(x) abs(squeeze(reshape(x(2,:,it0:it1),[],1)));
plot(plt(Ni00),plt(Ne00),'.');
title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
@@ -10,19 +10,19 @@ end
if 0
%% histograms
-t0 = 200; t1 = 300; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
+t0 = 2700; t1 = 2800; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
plt = @(x) squeeze(reshape(x(:,:,it0:it1),[],1));
-fig = figure;
+fig = figure('Position', [100, 100, 800, 500]);
subplot(221)
- hist(plt(ne00),100); hold on
+ hist(plt(ne00),50); hold on
title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
xlabel('$n_e$');
subplot(222)
- hist(plt(ni00),100); hold on
+ hist(plt(ni00),50); hold on
title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
xlabel('$n_i$');
subplot(223)
- hist(plt(phi),100); hold on
+ hist(plt(phi),50); hold on
title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
xlabel('$\phi$');
FIGNAME = ['histograms_ne_ni_phi'];
@@ -32,15 +32,15 @@ end
if 0
%% ni ne phi 3D phase space
fig = figure;
- t0 = 50; [~,it0] = min(abs(t0-Ts2D));
+ t0 = 2700; [~,it0] = min(abs(t0-Ts2D));
plt = @(x) squeeze(reshape(x(:,:,it0),[],1));
plot3(plt(ni00),plt(ne00),plt(phi),'.','MarkerSize',1.9);
title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
xlabel('$n_i$'); ylabel('$n_e$'); zlabel('$\phi$'); grid on;
end
-t0 = 0;
-skip_ = 1;
+t0 = 2400;
+skip_ = 2;
DELAY = 0.01*skip_;
FRAMES = floor(t0/(Ts2D(2)-Ts2D(1)))+1:skip_:numel(Ts2D);
if 0
diff --git a/wk/compute_collision_mat.m b/wk/compute_collision_mat.m
index 24841850538ee27f11590ceaad2ae8a400e70402..6a1f5e2842c7b67821fa6a4e4649d0e1468c5f05 100644
--- a/wk/compute_collision_mat.m
+++ b/wk/compute_collision_mat.m
@@ -2,7 +2,7 @@ addpath(genpath('../matlab')) % ... add
%% Grid configuration
N = 200; % Frequency gridpoints (Nkr = N/2)
L = 120; % Size of the squared frequency domain
-dk = 2*pi/L;
+dk = 2.*pi/L;
kmax = N/2*dk;
kr = dk*(0:N/2);
kz = dk*(0:N/2);
@@ -24,22 +24,22 @@ end
% end
%%
%% We compute only on a kperp grid with dk space from 0 to kperpmax
-kperp = unique([0:dk:(sqrt(2)*kmax),sqrt(2)*kmax]);
kperpmax = sqrt(2) * kmax;
+kperp = unique([0:dk:kperpmax,kperpmax]);
%%
n_ = 1;
for k_ = kperp
- disp(['----------Computing matrix ',num2str(n_),'/',num2str(Nperp),'----------'])
+ disp(['----------Computing matrix ',num2str(n_),'/',num2str(numel(kperp)),'----------'])
%% Script to run COSOlver in order to create needed collision matrices
COSOlver_path = '../../Documents/MoliSolver/COSOlver/';
- COSOLVER.pmaxe = 20;
- COSOLVER.jmaxe = 10;
- COSOLVER.pmaxi = 20;
- COSOLVER.jmaxi = 10;
+ COSOLVER.pmaxe = 10;
+ COSOLVER.jmaxe = 05;
+ COSOLVER.pmaxi = 10;
+ COSOLVER.jmaxi = 05;
COSOLVER.kperp = k_;
- COSOLVER.neFLR = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2))); % rule of thumb for sum truncation
- COSOLVER.niFLR = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2)));
+ COSOLVER.neFLR = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2)))+2; % rule of thumb for sum truncation
+ COSOLVER.niFLR = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2)))+2;
COSOLVER.idxT4max = 40;
COSOLVER.neFLRs = 0; % ... only for GK abel
@@ -63,14 +63,29 @@ for k_ = kperp
write_fort90_COSOlver
k_string = sprintf('%0.4f',k_);
- mat_file_name = ['../iCa/self_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
+ self_mat_file_name = ['../iCa/self_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
'_ESELF_',num2str(COSOLVER.co),'_ISELF_',num2str(COSOLVER.co),...
'_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
'_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
'_JE_12',...
'_NFLR_',num2str(COSOLVER.neFLR),...
'_kperp_',k_string,'.h5'];
- if (exist(mat_file_name,'file') > 0)
+ ei_mat_file_name = ['../iCa/ei_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
+ '_ETEST_',num2str(COSOLVER.co),'_EBACK_',num2str(COSOLVER.co),...
+ '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
+ '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
+ '_JE_12_tau_1.0000_mu_0.0233',...
+ '_NFLRe_',num2str(COSOLVER.neFLR),'_NFLRi_',num2str(COSOLVER.neFLR)...
+ '_kperp_',k_string,'.h5'];
+ ie_mat_file_name = ['../iCa/ie_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
+ '_ITEST_',num2str(COSOLVER.co),'_IBACK_',num2str(COSOLVER.co),...
+ '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
+ '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
+ '_JE_12_tau_1.0000_mu_0.0233',...
+ '_NFLRe_',num2str(COSOLVER.neFLR),'_NFLRi_',num2str(COSOLVER.neFLR)...
+ '_kperp_',k_string,'.h5'];
+% if (exist(self_mat_file_name,'file')>0 && exist(ei_mat_file_name,'file')>0 && exist(ie_mat_file_name,'file') > 0)
+ if (exist(ei_mat_file_name,'file')>0)%&& exist(ie_mat_file_name,'file') > 0)
disp(['Matrix available for kperp = ',k_string]);
else
cd ../../Documents/MoliSolver/COSOlver/
diff --git a/wk/continue_multiple_runs_marconi.m b/wk/continue_multiple_runs_marconi.m
new file mode 100644
index 0000000000000000000000000000000000000000..32202e1723f08691133a749766a3467fee9b2196
--- /dev/null
+++ b/wk/continue_multiple_runs_marconi.m
@@ -0,0 +1,74 @@
+%% Paste the list of continue_run calls
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.5_nu_1e+00_SGGK_CLOS_0_mu_2e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e+00_SGGK_CLOS_0_mu_2e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.7_nu_1e+00_SGGK_CLOS_0_mu_2e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.8_nu_1e+00_SGGK_CLOS_0_mu_2e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_5e-01_DGGK_CLOS_0_mu_2e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.7_nu_5e-01_DGGK_CLOS_0_mu_2e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.8_nu_5e-01_DGGK_CLOS_0_mu_2e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.9_nu_5e-01_DGGK_CLOS_0_mu_2e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.9_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt')
+
+%% Functions to modify preexisting fort.90 input file and launch on marconi
+function [] = continue_run(outfilename)
+ %% CLUSTER PARAMETERS
+ CLUSTER.PART = 'prod'; % dbg or prod
+ CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss
+ if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME = '00:30:00'; end;
+ CLUSTER.MEM = '64GB'; % Memory
+ CLUSTER.JNAME = 'HeLaZ';% Job name
+ NP_P = 2; % MPI processes along p
+ NP_KR = 24; % MPI processes along kr
+ % Compute processes distribution
+ Ntot = NP_P * NP_KR;
+ Nnodes = ceil(Ntot/48);
+ Nppn = Ntot/Nnodes;
+ CLUSTER.NODES = num2str(Nnodes); % MPI process along p
+ CLUSTER.NTPN = num2str(Nppn); % MPI process along kr
+ CLUSTER.CPUPT = '1'; % CPU per task
+ %%
+ RESDIR = ['../',outfilename(46:end-8),'/'];
+ BASIC.RESDIR = RESDIR;
+ FORT90 = [RESDIR,'fort.90'];
+
+ % Read txt into cell A
+ fid = fopen(FORT90,'r');
+ i = 1;
+ tline = fgetl(fid);
+ A{i} = tline;
+ while ischar(tline)
+ i = i+1;
+ tline = fgetl(fid);
+ A{i} = tline;
+ end
+ fclose(fid);
+
+ % Find previous job2load
+ if( numel(A{5}) == numel(' RESTART = .false.') )
+ A{5} = sprintf(' RESTART = .true.');
+ J2L = 0;
+ else
+ line = A{33};
+ line = line(end-2:end);
+ if(line(1) == '='); line = line(end); end;
+ J2L = str2num(line) + 1;
+ end
+ % Change job 2 load in fort.90
+ A{33} = [' job2load = ',num2str(J2L)];
+ disp(A{33})
+
+ % Rewrite fort.90
+ fid = fopen('fort.90', 'w');
+ for i = 1:numel(A)
+ if A{i+1} == -1
+ fprintf(fid,'%s', A{i});
+ break
+ else
+ fprintf(fid,'%s\n', A{i});
+ end
+ end
+ % Copy fort.90 into marconi
+ write_sbash_marconi
+ % Launch the job
+ system('ssh ahoffman@login.marconi.cineca.it sh HeLaZ/wk/setup_and_run.sh');
+end
\ No newline at end of file
diff --git a/wk/linear_study.m b/wk/linear_study.m
index 9734c39084ec0b5e81cd446fb8a7c0da2ae4ac8b..7b86263349924b29407c961dd2ad76987062828a 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -1,6 +1,6 @@
-for NU = [1.0 0.1 0.01]
-for ETAB = [0.5 0.6 0.7 0.8]
-for CO = [-3 -2 -1 0 1 2]
+for NU = [1.0]
+for ETAB = [0.5]
+for CO = [1,2]
%clear all;
addpath(genpath('../matlab')) % ... add
default_plots_options
@@ -14,11 +14,11 @@ TAU = 1.0; % e/i temperature ratio
% ETAB = 0.5;
ETAN = 1.0; % Density gradient
ETAT = 0.0; % Temperature gradient
-NU_HYP = 0.0; % Hyperdiffusivity coefficient
+NU_HYP = 1.0; % Hyperdiffusivity coefficient
LAMBDAD = 0.0;
NOISE0 = 1.0e-5;
%% GRID PARAMETERS
-N = 100; % Frequency gridpoints (Nkr = N/2)
+N = 200; % Frequency gridpoints (Nkr = N/2)
L = 120; % Size of the squared frequency domain
KREQ0 = 1; % put kr = 0
MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
@@ -62,13 +62,13 @@ MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
%% PARAMETER SCANS
if 1
%% Parameter scan over PJ
-PA = [2, 3, 4, 6, 8, 10];
-JA = [1, 2, 2, 3, 4, 5];
+% PA = [2, 3, 4, 6, 8, 10];
+% JA = [1, 2, 2, 3, 4, 5];
+PA = [4];
+JA = [2];
DTA= DT./sqrt(JA);
mup_ = MU_P;
muj_ = MU_J;
-% PA = [4];
-% JA = [1];
Nparam = numel(PA);
param_name = 'PJ';
gamma_Ni00 = zeros(Nparam,floor(N/2)+1);
@@ -88,15 +88,20 @@ for i = 1:Nparam
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 2 3; cd ../../../wk'])
% Load and process results
+ %%
load_results
tend = Ts2D(end); tstart = 0.4*tend;
+ [~,itstart] = min(abs(Ts2D-tstart));
+ [~,itend] = min(abs(Ts2D-tend));
for ikr = 1:N/2+1
- gamma_Ni00(i,ikr) = LinearFit_s(Ts2D,squeeze(abs(Ni00(ikr,1,:))),tstart,tend);
+ gamma_Ni00(i,ikr) = LinearFit_s(Ts2D(itstart:itend),squeeze(abs(Ni00(ikr,1,itstart:itend))));
Ni00_ST(i,ikr,1:numel(Ts2D)) = squeeze((Ni00(ikr,1,:)));
end
tend = Ts5D(end); tstart = 0.4*tend;
+ [~,itstart] = min(abs(Ts5D-tstart));
+ [~,itend] = min(abs(Ts5D-tend));
for ikr = 1:N/2+1
- gamma_Nipj(i,ikr) = LinearFit_s(Ts5D,squeeze(max(max(abs(Nipj(:,:,ikr,1,:)),[],1),[],2)),tstart,tend);
+ gamma_Nipj(i,ikr) = LinearFit_s(Ts5D(itstart:itend),squeeze(max(max(abs(Nipj(:,:,ikr,1,itstart:itend)),[],1),[],2)));
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));
diff --git a/wk/load_multiple_outputs_marconi.m b/wk/load_multiple_outputs_marconi.m
new file mode 100644
index 0000000000000000000000000000000000000000..955bbee4e4b9b4ff4da5913c3d49b45c01f86784
--- /dev/null
+++ b/wk/load_multiple_outputs_marconi.m
@@ -0,0 +1 @@
+%% Paste the list of simulation results to load
diff --git a/wk/local_run.m b/wk/local_run.m
index 6812cba4bdb741052d5298dc79aedbd0b1674983..b8dd706583ae8ff824ee9106691bed332fcdb9bb 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -4,37 +4,35 @@ addpath(genpath('../matlab')) % ... add
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 1.0; % Collision frequency
+NU = 0.1; % Collision frequency
ETAB = 0.6; % Magnetic gradient
ETAN = 1.0; % Density gradient
-NU_HYP = 1.0;
+NU_HYP = 0.2;
%% GRID PARAMETERS
N = 200; % Frequency gridpoints (Nkr = N/2)
-L = 80; % Size of the squared frequency domain
-PMAXE = 02; % Highest electron Hermite polynomial degree
-JMAXE = 01; % Highest '' Laguerre ''
-PMAXI = 02; % Highest ion Hermite polynomial degree
-JMAXI = 01; % Highest '' Laguerre ''
-MU_P = 0.0/PMAXI^2; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
-MU_J = 0.0/JMAXI^2; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
+L = 120; % Size of the squared frequency domain
+P = 2;
+J = 1;
+MU_P = 0.0/P^2; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
+MU_J = 0.0/J^2; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
%% TIME PARAMETERS
-TMAX = 2500; % Maximal time unit
-DT = 2e-2; % Time step
+TMAX = 500; % Maximal time unit
+DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for profiler
-SPS2D = 1/2; % Sampling per time unit for 2D arrays
-SPS5D = 1/4; % Sampling per time unit for 5D arrays
+SPS2D = 1; % Sampling per time unit for 2D arrays
+SPS5D = 1/50; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints/10
-RESTART = 1; % To restart from last checkpoint
-JOB2LOAD= 3;
+RESTART = 0; % To restart from last checkpoint
+JOB2LOAD= 0;
%% OPTIONS AND NAMING
% Collision operator
% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
-CO = 1;
-CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
+CO = 2;
+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_stability'; % Name of the simulation
-% SIMID = ['local_v2.5_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation
-% SIMID = sprintf(SIMID,NU);
+% SIMID = 'local_low_pj_analysis'; % Name of the simulation
+% SIMID = 'Sugama_no_ei_ie_coll'; % Name of the simulation
+SIMID = ['v2.5_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation
NON_LIN = 1; % activate non-linearity (is cancelled if KREQ0 = 1)
%% OUTPUTS
W_DOUBLE = 0;
@@ -46,6 +44,10 @@ W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% unused
+PMAXE = P; % Highest electron Hermite polynomial degree
+JMAXE = J; % Highest '' Laguerre ''
+PMAXI = P; % Highest ion Hermite polynomial degree
+JMAXI = J; % Highest '' Laguerre ''
KERN = 0; % Kernel model (0 : GK)
KR0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
KREQ0 = 0; % put kr = 0
@@ -61,4 +63,6 @@ ETAT = 0.0; % Temperature gradient
INIT_PHI= 1; % Start simulation with a noisy phi and moments
%% Setup and file management
setup
-system('rm fort.90');
\ No newline at end of file
+system('rm fort.90');
+outfile = [BASIC.RESDIR,'out.txt'];
+disp(outfile);
\ No newline at end of file
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 346e96e58546058a92f93e0de3803c4d1ffef279..868a20cbe9da870391c7cc6cf64331a5afb2472d 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -1,4 +1,4 @@
-%clear all;
+clear all;
addpath(genpath('../matlab')) % ... add
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
@@ -8,40 +8,40 @@ CLUSTER.PART = 'prod'; % dbg or prod
CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss
if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME = '00:30:00'; end;
CLUSTER.MEM = '64GB'; % Memory
-CLUSTER.JNAME = 'gamma_inf';% Job name
+CLUSTER.JNAME = 'HeLaZ';% Job name
NP_P = 2; % MPI processes along p
NP_KR = 24; % MPI processes along kr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.1; % Collision frequency
+NU = 1.0; % Collision frequency
ETAB = 0.6; % Magnetic gradient
NU_HYP = 1.0; % Hyperdiffusivity coefficient
+NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
+% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
+CO = -3;
%% GRID PARAMETERS
-N = 200; % Frequency gridpoints (Nkr = N/2)
-L = 80; % Size of the squared frequency domain
-P = 12; % Electron and Ion highest Hermite polynomial degree
-J = 06; % Electron and Ion highest Laguerre polynomial degree
-MU_P = 0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
-MU_J = 0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
+N = 200; % Frequency gridpoints (Nkr = N/2)
+L = 120; % Size of the squared frequency domain
+P = 10; % Electron and Ion highest Hermite polynomial degree
+J = 05; % Electron and Ion highest Laguerre polynomial degree
+MU_P = 0.0/(P/2)^4;% Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
+MU_J = 0.0/(J/2)^2;% Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
%% TIME PARAMETERS
-TMAX = 4000; % Maximal time unit
+TMAX = 5000; % Maximal time unit
DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1/2; % 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
RESTART = 1; % To restart from last checkpoint
-JOB2LOAD= 4;
-%% OPTIONS
-SIMID = ['HeLaZ_v2.5_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation
-SIMID = sprintf(SIMID,NU);
-% SIMID = 'test_marconi_sugama'; % Name of the simulation
+JOB2LOAD= 0;
+%% Naming
+% SIMID = 'test'; % Name of the simulation
+SIMID = ['v2.5_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation
PREFIX =[];
% PREFIX = sprintf('%d_%d_',NP_P, NP_KR);
-% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
-CO = 1;
+%% Options
CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
-NL_CLOS = 1; % nonlinear closure model (0: =0 nmax = jmax, 1: nmax = jmax-j, >1 : nmax = NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
INIT_PHI= 1; % Start simulation with a noisy phi and moments
%% OUTPUTS
@@ -64,6 +64,7 @@ JMAXE = J; % Highest '' Laguerre ''
PMAXI = P; % Highest ion Hermite polynomial degree
JMAXI = J; % Highest '' Laguerre ''
kmax = N*pi/L;% Highest fourier mode
+% kmax = 2/3*N*pi/L;% Highest fourier mode with AA
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
NOISE0 = 1.0e-5;
diff --git a/wk/open_figure_script.m b/wk/open_figure_script.m
new file mode 100644
index 0000000000000000000000000000000000000000..2b82c3a0df89b2c8f6be6343e5d5c7615955a2a4
--- /dev/null
+++ b/wk/open_figure_script.m
@@ -0,0 +1,24 @@
+simname_ = '';
+
+%%
+% simname_ = 'v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.7_nu_1e+00_SGGK_CLOS_0_mu_2e-02';
+% simname_ = 'v2.5_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.8_nu_1e+00_SGGK_CLOS_0_mu_2e-02';
+
+% simname_ = 'v2.5_P_2_J_1/200x100_L_120_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02';
+% simname_ = 'v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02';
+simname_ = 'v2.5_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_5e-01_DGGK_CLOS_0_mu_2e-02';
+
+
+%%
+% figname_ = '/fig/ZF_transport_drphi_';
+% figname_ = '/fig/space_time_';
+figname_ = '/fig/phi_shear_phase_space_';
+
+path_ = '../results/';
+
+[~,idx] = max(simname_=='x');
+
+params_ = simname_(idx-3:end);
+
+
+openfig([path_,simname_,figname_,params_,'.fig']);
\ No newline at end of file
diff --git a/wk/spectral_analysis.m b/wk/spectral_analysis.m
new file mode 100644
index 0000000000000000000000000000000000000000..afb9f7177289308e71a7e05ae775449d7cd4d4ea
--- /dev/null
+++ b/wk/spectral_analysis.m
@@ -0,0 +1,67 @@
+
+%%
+if 0
+%% Hermite energy spectra
+% tf = Ts2D(end-3);
+skip = 1;
+fig = figure; FIGNAME = ['hermite_spectrum_',PARAMS];set(gcf, 'Position', [100, 100, 1800, 600]);
+plt = @(x) squeeze(x);
+N10 = floor(Nji/10);
+Nrow = (1+N10);
+Ncol = ceil(Nji/Nrow);
+for ij = 1:Nji
+ subplot(Nrow,Ncol,ij)
+ for it5 = 1:skip:Ns5D
+ alpha = it5*1.0/Ns5D;
+ loglog(Pi(1:2:end),plt(epsilon_i_pj(1:2:end,ij,it5)),...
+ 'color',(1-alpha)*[0.8500, 0.3250, 0.0980]+alpha*[0, 0.4470, 0.7410],...
+ 'DisplayName',['t=',num2str(Ts5D(it5))]); hold on;
+ end
+ grid on; ylim([1e0,1e10]);
+ xlabel('$p$');
+ TITLE = ['$\sum_{kr,kz} |N_i^{p',num2str(Ji(ij)),'}|^2$']; title(TITLE);
+end
+save_figure
+end
+%%
+if 0
+%% Laguerre energy spectra
+% tf = Ts2D(end-3);
+skip = 1;
+fig = figure; FIGNAME = ['laguerre_spectrum_',PARAMS];set(gcf, 'Position', [100, 100, 1800, 600]);
+plt = @(x) squeeze(x);
+N10 = floor(Npi/10);
+Nrow = ceil((1+N10)/2);
+Ncol = ceil(Npi/Nrow/2);
+for ip = 1:2:Npi
+ subplot(Nrow,Ncol,ip/2+1)
+ for it5 = 1:skip:Ns5D
+ alpha = it5*1.0/Ns5D;
+ loglog(Ji,plt(epsilon_i_pj(ip,:,it5)),...
+ 'color',(1-alpha)*[0.8500, 0.3250, 0.0980]+alpha*[0, 0.4470, 0.7410],...
+ 'DisplayName',['t=',num2str(Ts5D(it5))]); hold on;
+ grid on;
+ xlabel('$j$'); ylim([1e-20,1e10]);
+ TITLE = ['$\sum_{kr,kz} |N_i^{',num2str(Pi(ip)),'j}|^2$']; title(TITLE);
+ end
+end
+save_figure
+end
+
+%%
+no_AA = (2:floor(2*Nkr/3));
+tKHI = 100;
+[~,itKHI] = min(abs(Ts2D-tKHI));
+after_KHI = (itKHI:Ns2D);
+if 0
+%% Phi frequency space time diagram at kz=kz(ikz)
+kz_ = 0.0;
+[~,ikz] = min(abs(kz-kz_));
+fig = figure; FIGNAME = ['phi_freq_diag_',PARAMS];set(gcf, 'Position', [100, 100, 500, 400]);
+ [TY,TX] = meshgrid(Ts2D(after_KHI),kr(no_AA));
+ pclr = pcolor(TX,TY,(squeeze(abs(PHI(no_AA,ikz,(after_KHI)))))); set(pclr, 'edgecolor','none'); colorbar;
+ caxis([0,10000]); colormap hot
+ ylabel('$t c_s/R$'), xlabel('$0<k_r<2/3 k_r^{\max}$')
+ legend(['$|\tilde\phi(k_z=',num2str(kz_),')|$'])
+ title('Spectrogram of $\phi$')
+end
\ No newline at end of file
diff --git a/wk/transport_results_2_5.m b/wk/transport_results_2_5.m
new file mode 100644
index 0000000000000000000000000000000000000000..c47aca0dedb0a2a246cee4751a70f0eeb5084db8
--- /dev/null
+++ b/wk/transport_results_2_5.m
@@ -0,0 +1,103 @@
+%% nuDGGK = 1.0
+figure
+% 6,3 nuDGGK=1.0
+eta = [0.5,0.6,0.7,0.8];
+gamma = [32.6443 3.6895 0.3744 0.056];
+std_ = [6.1529 0.7986 0.0493 0.0134];
+errorbar(eta,gamma,std_,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=1.0$',...
+ 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(1,:)); hold on;
+
+% 10,5 nuDGGK = 1.0
+eta = [0.6,0.7,0.8];
+gamma = [3.5546 0.3917 0.0421];
+std_ = [0.5846 0.0486 0.0203];
+errorbar(eta,gamma,std_,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=1.0$',...
+ 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(2,:)); hold on;
+
+% 12,6 nuDGGK = 1.0
+eta = [0.6];
+gamma = [4.064];
+std_ = [0.7964];
+errorbar(eta,gamma,std_,'-','DisplayName','$P,J=12,06$, $\nu_{DGGK}=1.0$',...
+ 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(3,:)); hold on;
+
+
+% Mix len Ricci 2006
+semilogy([0.5 1.0],[10 1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
+set(gca, 'YScale', 'log'); grid on; legend('show')
+xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$')
+
+%% nuDGGK = 0.5
+figure
+% 6,3 nuDGGK=0.5
+eta = [0.5,0.6,0.7,0.8];
+gamma = [20.511 2.6292 0.2353 0.057];
+std_ = [3.67 1.13 0.055 0.004];
+errorbar(eta,gamma,std_,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.5$',...
+ 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(4,:)); hold on;
+
+% 10,5 nuDGGK=0.5
+eta = [0.6,0.7,0.8];
+gamma = [2.4949 0.2578 0.058];
+std_ = [0.8696 0.1191 0.011];
+errorbar(eta,gamma,std_,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.5$',...
+ 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(5,:)); hold on;
+
+% Mix len Ricci 2006
+semilogy([0.5 1.0],[10 1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
+set(gca, 'YScale', 'log'); grid on; legend('show')
+xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$')
+
+%% nuDGGK = 0.1
+figure
+% 6,3
+eta = [0.6 0.7 0.8];
+gamma = [0.24321 0.085 0.0367];
+std_ = [0.295 0.05 0.0023];
+bmax = [1.0 0.21 0.0367];
+bmin = [0.02 0.04 0.0367];
+% errorbar(eta,gamma,std_,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.1$',...
+% 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(6,:)); hold on;
+plot(eta,gamma,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
+ 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(6,:)); hold on;
+plot(eta,bmax,'-.','DisplayName',' ','Color',line_colors(6,:)); hold on;
+plot(eta,bmin,'-.','DisplayName',' ','Color',line_colors(6,:)); hold on;
+
+% 10,5
+eta = [0.5 0.6 0.7 0.8 0.9];
+gamma = [12.2 0.2133 0.088 0.04253 0.02435];
+std_ = [4.7 0.21328 0.065 0.009 0.00386];
+bmax = [12.2 0.8 0.25 0.04253 0.02435];
+bmin = [12.2 0.02 0.03 0.04253 0.02435];
+% errorbar(eta,gamma,std_,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
+% 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
+plot(eta,gamma,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
+ 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
+plot(eta,bmax,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
+plot(eta,bmin,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
+% Mix len Ricci 2006
+semilogy([0.5 1.0],[10 1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
+set(gca, 'YScale', 'log'); grid on; legend('show')
+xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$')
+
+%% nuSGGK = 1.0
+figure
+% 6,3 nuDGGK=1.0
+eta = [0.5 0.6,0.7,0.8];
+gamma = [2.3 0.2215 0.0133 0.0032];
+std_ = [3.1 0.22 0.0019 0.0006];
+errorbar(eta,gamma,std_,'-','DisplayName','$P,J=06,03$, $\nu_{SGGK}=1.0$',...
+ 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(8,:)); hold on;
+
+
+% 10,5 nuDGGK = 1.0
+eta = [0.5 0.6,0.7,0.8];
+gamma = [10 0.319 0.009 0.0026];
+std_ = [1.34 0.228 0.001 0.001];
+% errorbar(eta,gamma,err,'-.','DisplayName','$P,J=10,5$, $\nu_{SGGK}=1.0$',...
+% 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(1,:)); hold on;
+
+% Mix len Ricci 2006
+semilogy([0.5 1.0],[10 1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
+set(gca, 'YScale', 'log'); grid on; legend('show')
+xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$')