From a7017e860f275099715fc3b67dba1670d3ddce02 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Tue, 25 Apr 2023 11:45:11 +0200
Subject: [PATCH] save scripts
---
matlab/assemble_cosolver_matrices.m | 71 ++++++++
matlab/compile_results.m | 4 +-
matlab/compute/mode_growth_meter.m | 10 +-
matlab/extract_fig_data.m | 3 +-
matlab/load/compile_results_low_mem.m | 4 +-
matlab/load/load_write_mat.m | 15 ++
matlab/load/quick_gene_load.m | 5 +-
matlab/plot/create_film.m | 2 +-
matlab/plot/plot_cosol_mat.m | 99 +++++++----
.../plot_radial_transport_and_spacetime.m | 7 +-
matlab/plot/show_moments_spectrum.m | 10 +-
matlab/plot/top_title.m | 9 +
matlab/setup.m | 23 +--
matlab/write_fort90.m | 1 -
.../old => }/CBC_kT_nu_scan.m | 114 ++++++------
wk/analysis_gene.m | 29 +++-
.../GX_2022_fig_1_1.txt | 11 ++
.../GX_2022_fig_2_1.txt | 23 +++
.../GX_2022_fig_4_1.txt | 93 ++++++++++
.../GX_2022_fig_4_2.txt | 5 +
wk/benchmark_and_scan_scripts/old/fort_00.90 | 92 ++++++++++
wk/benchmark_and_scan_scripts/plot_GX_data.m | 33 ++++
wk/fast_analysis.m | 164 ++++++++----------
wk/fast_heat_flux_analysis.m | 43 +++++
wk/heat_flux_convergence_analysis.m | 156 +++++++++++++++++
wk/kTcrit_convergence_analysis.m | 146 ++++++++++++++++
wk/lin_ITG_salpha.m | 162 +++++++++++++++++
wk/lin_Ivanov.m | 5 +-
wk/load_metadata_scan.m | 17 +-
wk/local_run.m | 38 ++--
wk/multiple_kTscan_p2_analysis.m | 113 +++++++++---
wk/old scripts/header_dp.m | 67 +++++++
wk/scan_lin_ITG.m | 163 +++++++++++++++++
33 files changed, 1474 insertions(+), 263 deletions(-)
create mode 100644 matlab/assemble_cosolver_matrices.m
create mode 100644 matlab/load/load_write_mat.m
create mode 100644 matlab/plot/top_title.m
rename wk/{benchmark_and_scan_scripts/old => }/CBC_kT_nu_scan.m (64%)
create mode 100644 wk/benchmark_and_scan_scripts/GX_2022_fig_1_1.txt
create mode 100644 wk/benchmark_and_scan_scripts/GX_2022_fig_2_1.txt
create mode 100644 wk/benchmark_and_scan_scripts/GX_2022_fig_4_1.txt
create mode 100644 wk/benchmark_and_scan_scripts/GX_2022_fig_4_2.txt
create mode 100644 wk/benchmark_and_scan_scripts/old/fort_00.90
create mode 100644 wk/benchmark_and_scan_scripts/plot_GX_data.m
create mode 100644 wk/fast_heat_flux_analysis.m
create mode 100644 wk/heat_flux_convergence_analysis.m
create mode 100644 wk/kTcrit_convergence_analysis.m
create mode 100644 wk/lin_ITG_salpha.m
create mode 100644 wk/old scripts/header_dp.m
create mode 100644 wk/scan_lin_ITG.m
diff --git a/matlab/assemble_cosolver_matrices.m b/matlab/assemble_cosolver_matrices.m
new file mode 100644
index 00000000..bf71b20b
--- /dev/null
+++ b/matlab/assemble_cosolver_matrices.m
@@ -0,0 +1,71 @@
+codir = '/home/ahoffman/cosolver/';
+matname= 'gk_sugama_tau1e-3_P4_J2_dk_5e-2_km_5.0_NFLR_5';
+matdir = [codir,'SGGK_tau1e-3_P4_J2_dk_0.05_kpm_5.0/matrices/'];
+kperp = load([matdir,'kperp.in']);
+Nkp = numel(kperp);
+outdir = '/home/ahoffman/gyacomo/iCa/';
+outfile= [outdir,matname,'.h5'];
+
+if(exist(outfile)>0)
+ system(['rm ',outdir,matname,'.h5']);
+end
+
+Nmax = numel(kperp);
+for n_=1:Nmax
+ n_string = sprintf('%5.5d',n_); disp(n_string);
+ filename = ['ei.',n_string,'.h5'];
+ if(exist([matdir,filename])>0)
+ % Load matrices and write them in one h5 file
+
+ olddname = '/Ceipj/CeipjF'; infile = [matdir,'ei.',n_string,'.h5'];
+ newdname = ['/',sprintf('%5.5d',n_-1),olddname];
+ load_write_mat(infile,olddname,outfile,newdname);
+
+ olddname = '/Ceipj/CeipjT'; infile = [matdir,'ei.',n_string,'.h5'];
+ newdname = ['/',sprintf('%5.5d',n_-1),olddname];
+ load_write_mat(infile,olddname,outfile,newdname);
+
+ olddname = '/Ciepj/CiepjT'; infile = [matdir,'ie.',n_string,'.h5'];
+ newdname = ['/',sprintf('%5.5d',n_-1),olddname];
+ load_write_mat(infile,olddname,outfile,newdname);
+
+ olddname = '/Ciepj/CiepjF'; infile = [matdir,'ie.',n_string,'.h5'];
+ newdname = ['/',sprintf('%5.5d',n_-1),olddname];
+ load_write_mat(infile,olddname,outfile,newdname);
+
+ olddname = '/Caapj/Ceepj'; infile = [matdir,'self.',n_string,'.h5'];
+ newdname = ['/',sprintf('%5.5d',n_-1),olddname];
+ mat_ee = load_write_mat(infile,olddname,outfile,newdname);
+
+ olddname = '/Caapj/Ciipj'; infile = [matdir,'self.',n_string,'.h5'];
+ newdname = ['/',sprintf('%5.5d',n_-1),olddname];
+ mat_ii = load_write_mat(infile,olddname,outfile,newdname);
+
+ % to verify symmetry
+ verif = max(abs(imag(eig(mat_ee)))) + max(abs(imag(eig(mat_ii))));
+ if(verif>0) disp(['Warning, non sym matrix at ', n_string]); end;
+ % to verify negative EV
+ verif = max(real(eig(mat_ee))) + max(real(eig(mat_ii)));
+ if(verif>0) disp(['Warning, positive EV at ', n_string]); end;
+ else
+ break
+ end
+end
+olddname = '/Ceipj/CeipjF'; infile = [matdir,'ei.',n_string,'.h5'];
+Pmaxe = h5readatt(infile,olddname,'Pmaxe');
+Jmaxe = h5readatt(infile,olddname,'Jmaxe');
+Pmaxi = h5readatt(infile,olddname,'Pmaxi');
+Jmaxi = h5readatt(infile,olddname,'Jmaxi');
+dims_e = [0 0];
+dims_e(1)= Pmaxe; dims_e(2) = Jmaxe;
+dims_i = [Pmaxi; Jmaxi];
+h5create(outfile,'/dims_e',numel(dims_e));
+h5write (outfile,'/dims_e',dims_e);
+h5create(outfile,'/dims_i',numel(dims_i));
+h5write (outfile,'/dims_i',dims_i);
+%
+h5create(outfile,'/coordkperp',numel(kperp(1:Nmax)));
+h5write (outfile,'/coordkperp',kperp(1:Nmax)');
+fid = H5F.open(outfile);
+
+H5F.close(fid);
\ No newline at end of file
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 27c14463..e744f989 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -60,8 +60,8 @@ while(CONTINUE)
%% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(sprintf('Loading ID %.2d (%s)',JOBNUM,filename));
% Loading from output file
- CPUTIME = h5readatt(filename,'/data/input','cpu_time');
- DT_SIM = h5readatt(filename,'/data/input','dt');
+ % CPUTIME = h5readatt(filename,'/data/input','cpu_time');
+ % 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');
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index 30ab59ff..907aa85e 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -10,11 +10,11 @@ if numel(size(DATA.PHI)) == 3
else
switch OPTIONS.iz
case 'avg'
- field = squeeze(mean(DATA.PHI,3));
+ field = reshape(mean(DATA.PHI,3),DATA.grids.Nky,DATA.grids.Nkx,numel(DATA.Ts3D));
zstrcomp = 'z-avg';
otherwise
- field = squeeze(DATA.PHI(:,:,OPTIONS.iz,:));
- zstrcomp = ['z=',DATA.z(OPTIONS.iz)];
+ field = reshape(DATA.PHI(:,:,OPTIONS.iz,:),DATA.grids.Nky,DATA.grids.Nkx,numel(DATA.Ts3D));
+ zstrcomp = ['z=',num2str(DATA.grids.z(OPTIONS.iz))];
end
end
@@ -30,7 +30,7 @@ TW = [OPTIONS.KY_TW; OPTIONS.KX_TW];
[~,ikzf] = max(mean(squeeze(abs(field(1,:,FRAMES))),2));
-FIGURE.fig = figure; set(gcf, 'Position', [100 100 1200 700])
+FIGURE.fig = figure; %set(gcf, 'Position', [100 100 1200 700])
FIGURE.FIGNAME = 'mode_growth_meter';
for i = 1:2
MODES_SELECTOR = i; %(1:kx=0; 2:ky=0)
@@ -185,5 +185,5 @@ if d
title('Growth rates')
end
-suptitle(DATA.paramshort)
+top_title(DATA.paramshort)
end
\ No newline at end of file
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index 825aff05..a9689411 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -19,7 +19,8 @@ end
% n1 = numel(X_);
[~,n0] = min(abs(X_-tw(1)));
[~,n1] = min(abs(X_-tw(2)));
-mvm = @(x) movmean(x,1);
+mvm = @(x) x;
+% mvm = @(x) movmean(x,1);
shift = 0;%X_(n0);
% shift = 0;
skip = 1;
diff --git a/matlab/load/compile_results_low_mem.m b/matlab/load/compile_results_low_mem.m
index dcc9ad2b..54ee00f4 100644
--- a/matlab/load/compile_results_low_mem.m
+++ b/matlab/load/compile_results_low_mem.m
@@ -34,7 +34,7 @@ while(CONTINUE)
end
if openable
%% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- disp(sprintf('Loading ID %.2d (%s)',JOBNUM,filename));
+ fprintf('Loading ID %.2d (%s)\n',JOBNUM,filename);
% Loading from output file
CPUTIME = h5readatt(filename,'/data/input','cpu_time');
DT_SIM = h5readatt(filename,'/data/input/basic','dt');
@@ -148,7 +148,7 @@ else
if(numel(Parray)>1)
DATA.grids.deltap = Parray(2)-Parray(1);
else
- DATA.grids.deltap = 1
+ DATA.grids.deltap = 1;
end
DATA.dir = DIRECTORY;
DATA.localdir = DIRECTORY;
diff --git a/matlab/load/load_write_mat.m b/matlab/load/load_write_mat.m
new file mode 100644
index 00000000..3315b172
--- /dev/null
+++ b/matlab/load/load_write_mat.m
@@ -0,0 +1,15 @@
+function mat_ = load_write_mat(infile,olddname,outfile,newdname)
+ mat_ = h5read (infile,olddname);
+
+ h5create (outfile,newdname,size(mat_));
+ h5write (outfile,newdname,mat_);
+
+ if(~strcmp(olddname(end-3:end-2),'ii'))
+ neFLR = h5readatt(infile,olddname,'neFLR');
+ h5writeatt(outfile,newdname,'neFLR',neFLR);
+ end
+ if(~strcmp(olddname(end-3:end-2),'ee'))
+ niFLR = h5readatt(infile,olddname,'niFLR');
+ h5writeatt(outfile,newdname,'niFLR',niFLR);
+ end
+end
diff --git a/matlab/load/quick_gene_load.m b/matlab/load/quick_gene_load.m
index 3ca0baa1..ef4d4c2e 100644
--- a/matlab/load/quick_gene_load.m
+++ b/matlab/load/quick_gene_load.m
@@ -77,6 +77,9 @@ path = '/home/ahoffman/gene/linear_CBC_results/';
% fname = 'CBC_miller_nz_24_nv_scan_nw_16_kine.txt';
%---------- CBC
% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
+% fname = 'CBC_salpha_nx_8_nz_18_nv_18_nw_8_adiabe.txt';
+% fname = 'CBC_salpha_nx_8_nz_18_nv_12_nw_8_adiabe.txt';
+fname = 'CBC_salpha_nx_8_nz_24_nv_60_nw_30_adiabe.txt';
% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_kine.txt';
% fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_adiabe.txt';
% fname = 'CBC_miller_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
@@ -85,7 +88,7 @@ path = '/home/ahoffman/gene/linear_CBC_results/';
%----------Convergence nv CBC
% fname = 'CBC_ky_0.3_nv_scan_8x1x24_nw_8_Lv_3_Lw_6.txt';
% fname = 'CBC_ky_0.3_nv_scan_8x1x24_nw_16_Lv_3_Lw_6.txt';
-fname = 'CBC_ky_0.3_nv_scan_8x1x24_nw_24_Lv_3_Lw_6.txt';
+% fname = 'CBC_ky_0.3_nv_scan_8x1x24_nw_24_Lv_3_Lw_6.txt';
data_ = load([path,fname]);
diff --git a/matlab/plot/create_film.m b/matlab/plot/create_film.m
index 58b1baef..e727e49c 100644
--- a/matlab/plot/create_film.m
+++ b/matlab/plot/create_film.m
@@ -36,7 +36,7 @@ if hmax == hmin
disp('Warning : h = hmin = hmax = const')
else
% Setup figure frame
-fig = figure('Color','white','Position', toplot.DIMENSIONS.*[1 1 1.2 1]);
+fig = figure('Color','white','Position', toplot.DIMENSIONS.*[0 0 1.2 1]);
if ~strcmp(OPTIONS.PLAN,'sx')
pcolor(X,Y,FIELD(:,:,1)); % to set up
if BWR
diff --git a/matlab/plot/plot_cosol_mat.m b/matlab/plot/plot_cosol_mat.m
index 686686f8..65e04b31 100644
--- a/matlab/plot/plot_cosol_mat.m
+++ b/matlab/plot/plot_cosol_mat.m
@@ -97,30 +97,46 @@ subplot(224)
if 0
%%
mfns = {...
- '/home/ahoffman/gyacomo/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
-% '/home/ahoffman/gyacomo/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
-% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_4.h5',...
-% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12.h5',...
-% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5',...
-% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5',...
- '/home/ahoffman/gyacomo/iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5',...
-% '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
- '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
-% '/home/ahoffman/gyacomo/iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5',...
-% '/home/ahoffman/gyacomo/iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5',...
+ '/home/ahoffman/gyacomo/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5';...
+% '/home/ahoffman/gyacomo/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5';...
+% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_4.h5';...
+% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12.h5';...
+% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5';...
+% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5';...
+ '/home/ahoffman/gyacomo/iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5';...
+% '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5';...
+ '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';...
+ '/home/ahoffman/gyacomo/iCa/gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5';...
+ '/home/ahoffman/gyacomo/iCa/gk_landau_P11_J7_dk_5e-2_km_2.0_NFLR_16.h5';...
+ '/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-3_P11_J7_dk_5e-2_km_5.0_NFLR_12.h5';...
+ '/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-3_P4_J2_dk_5e-2_km_5.0_NFLR_5.h5';...
+% '/home/ahoffman/gyacomo/iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5';...
+% '/home/ahoffman/gyacomo/iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5';...
};
CONAME_A = {...
- 'SG 20 10',...
-% 'PA 20 10',...
-% 'FC 10 5 NFLR 4',...
-% 'FC 10 5 NFLR 12',...
-% 'FC 10 5 NFLR 12 k<2', ...
- 'LD 6 3 NFLR 20', ...
-% 'FC 4 2 NFLR 6',...
- 'FC 4 2 NFLR 12', ...
-% 'Hacked SG A',...
-% 'Hacked SG B',...
+ 'SG 20 10';...
+% 'PA 20 10';...
+% 'FC 10 5 NFLR 4';...
+% 'FC 10 5 NFLR 12';...
+% 'FC 10 5 NFLR 12 k<2';, ...
+ 'LD 6 3 NFLR 20'; ...
+% 'FC 4 2 NFLR 6';...
+ 'FC 4 2 NFLR 12'; ...
+ 'LD 10 5 NFLR 12'; ...
+ 'LD 11 7 NFLR 16'; ...
+ 'SG 11 7 NFLR 12, tau 1e-3'; ...
+ 'SG 4 2 NFLR 5, tau 1e-3'; ...
+% 'Hacked SG A';...
+% 'Hacked SG B';...
};
+TAU_A = [1;...
+ 1;...
+ 1;...
+ 1;...
+ 1;...
+ 1e-3;...
+ 1e-3;...
+ ];
figure
for j_ = 1:numel(mfns)
mat_file_name = mfns{j_}; disp(mat_file_name);
@@ -135,9 +151,9 @@ for j_ = 1:numel(mfns)
wmax(idx_+1) = max(abs(imag(eig(MAT))));
end
subplot(121)
- plot(kp_a,gmax,'DisplayName',CONAME_A{j_}); hold on;
+ plot(kp_a,gmax*TAU_A(j_),'DisplayName',CONAME_A{j_}); hold on;
subplot(122)
- plot(kp_a,wmax,'DisplayName',CONAME_A{j_}); hold on;
+ plot(kp_a,wmax*TAU_A(j_),'DisplayName',CONAME_A{j_}); hold on;
end
subplot(121)
legend('show'); grid on;
@@ -152,13 +168,34 @@ end
if 0
%%
kperp= 1.5;
-mfns = {'/home/ahoffman/gyacomo/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
- '/home/ahoffman/gyacomo/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
- '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
- '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
- '/home/ahoffman/gyacomo/iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5',...
+mfns = {...
+ '/home/ahoffman/gyacomo/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5';...
+ '/home/ahoffman/gyacomo/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5';...
+ '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5';...
+ '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';...
+ '/home/ahoffman/gyacomo/iCa/gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5';...
+ '/home/ahoffman/gyacomo/iCa/gk_landau_P11_J7_dk_5e-2_km_2.0_NFLR_16.h5';...
+ '/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-3_P11_J7_dk_5e-2_km_5.0_NFLR_12.h5';...
+ '/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-3_P4_J2_dk_5e-2_km_5.0_NFLR_5.h5';...
};
-CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'LD 6 3 NFLR 12'};
+CONAME_A = {'SG 20 10';...
+ 'PA 20 10';...
+ 'FC 4 2 NFLR 6';...
+ 'FC 4 2 NFLR 12';...
+ 'LD 10 5 NFLR 12';...
+ 'LD 11 7 NFLR 16';...
+ 'SG 11 7 NFLR 12, tau 1e-3'; ...
+ 'SG 4 2 NFLR 5, tau 1e-3'; ...
+ };
+TAU_A = [1;...
+ 1;...
+ 1;...
+ 1;...
+ 1;...
+ 1;...
+ 1e-3;...
+ 1e-3;...
+ ];
grow = {};
puls = {};
for j_ = 1:numel(mfns)
@@ -167,8 +204,8 @@ for j_ = 1:numel(mfns)
[~,idx_] = min(abs(kp_a-kperp));
matidx = sprintf('%5.5i',idx_);
MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
- grow{j_} = real(eig(MAT));
- puls{j_} = imag(eig(MAT));
+ grow{j_} = real(eig(MAT))*TAU_A(j_);
+ puls{j_} = imag(eig(MAT))*TAU_A(j_);
end
figure
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 84dacb44..4e2617cd 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -28,10 +28,10 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
FIGURE.ax1 = subplot(2,1,1,'parent',FIGURE.fig);
for ia = 1:DATA.inputs.Na
plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI(ia,:)*SCALE),'--',...
- 'color',clr_((DATA.grids.Np-1)/2+(ia-1),:),...
+ 'color',clr_(max(1,(DATA.grids.Np-1)/2+(ia-1)),:),...
'DisplayName',['$\Gamma_x$ ',DATA.paramshort]); hold on;
plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X(ia,:)*SCALE),'-',...
- 'color',clr_((DATA.grids.Np-1)/2+(ia-1),:),...
+ 'color',clr_(max(1,(DATA.grids.Np-1)/2+(ia-1)),:),...
'DisplayName',['$Q_x$ ',DATA.paramshort]); hold on;
ylabel('Transport')
if(~isnan(Qx_infty_avg))
@@ -84,6 +84,5 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
subplot(311)
plot(DATA.Ts3D,squeeze(mean(plt(f2plot),1)));
end
- suptitle(DATA.paramshort)
-
+ top_title(DATA.paramshort)
end
\ No newline at end of file
diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m
index 3dcba78b..b77806a7 100644
--- a/matlab/plot/show_moments_spectrum.m
+++ b/matlab/plot/show_moments_spectrum.m
@@ -5,18 +5,16 @@ Pa = DATA.grids.Parray;
Ja = DATA.grids.Jarray;
Time_ = DATA.Ts3D;
FIGURE.fig = figure; FIGURE.FIGNAME = ['mom_spectrum_',DATA.params_string];
-set(gcf, 'Position', [100 50 1000 400])
+% set(gcf, 'Position', [100 50 1000 400])
if OPTIONS.ST == 0
OPTIONS.LOGSCALE = 0;
end
if OPTIONS.LOGSCALE
logname = 'log';
-% compress = @(x,ia) (log(sum(abs(x(ia,:,:,:,:)),4)));
- compress = @(x,ia) log(sum(abs(squeeze(x(:,:,:,:))),3));
+ compress = @(x,ia) log(sum(abs((x(:,:,:,:,:))),4));
else
logname = '';
-% compress = @(x,ia) (sum(abs(x(ia,:,:,:,:)),4));
- compress = @(x,ia) sum(abs(squeeze(x(:,:,:,:))),3);
+ compress = @(x,ia) sum(abs((x(:,:,:,:,:))),4);
end
for ia = 1:DATA.inputs.Na
Napjz = compress(DATA.Napjz,ia);
@@ -99,7 +97,7 @@ for ia = 1:DATA.inputs.Na
end
title(plotname)
end
-suptitle(DATA.paramshort)
+top_title(DATA.paramshort)
end
diff --git a/matlab/plot/top_title.m b/matlab/plot/top_title.m
new file mode 100644
index 00000000..5e59141d
--- /dev/null
+++ b/matlab/plot/top_title.m
@@ -0,0 +1,9 @@
+function [] = top_title(title)
+% Home made function to use suptitle or sgtitle
+try % valid for matlab 2016
+ suptitle(title)
+catch % valid for matlab 2023
+ sgtitle(title)
+end
+
+end
\ No newline at end of file
diff --git a/matlab/setup.m b/matlab/setup.m
index 6256ca69..3236997d 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -61,25 +61,28 @@ CLOSURE.nmax = NMAX;
COLL.collision_model = ['''',CO,''''];
if (GKCO); COLL.GK_CO = '.true.'; else; COLL.GK_CO = '.false.';end;
if (ABCO); COLL.INTERSPECIES = '.true.'; else; COLL.INTERSPECIES = '.false.';end;
-COLL.mat_file = '''null''';
+COLL.mat_file = 'null';
switch CO
case 'SG'
- COLL.mat_file = '''../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5''';
-% COLL.mat_file = '''../../../iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5''';
-% COLL.mat_file = '''../../../iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5''';
+ COLL.mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5';
+% COLL.mat_file = 'gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5';
+% COLL.mat_file = 'gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5';
case 'LR'
- COLL.mat_file = '''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5''';
+ COLL.mat_file = 'gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5';
case 'LD'
- COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5''';
-% COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5''';
-% COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5''';
-% COLL.mat_file = '''../../../iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5''';
+ COLL.mat_file = 'gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5';
+ % COLL.mat_file = 'gk_landau_P11_J7_dk_5e-2_km_2.0_NFLR_16.h5';
+ % COLL.mat_file = 'gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';
+% COLL.mat_file = 'LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5';
+% COLL.mat_file = 'LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5';
+% COLL.mat_file = 'LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5';
end
+COLL.mat_file = [gyacomodir,'iCa/',COLL.mat_file];
+COLL.mat_file = ['''',COLL.mat_file,''''];
COLL.coll_kcut = COLL_KCUT;
% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = ['''',NUMERICAL_SCHEME,''''];
INITIAL.INIT_OPT = ['''',INIT_OPT,''''];
-INITIAL.ACT_ON_MODES = ['''',ACT_ON_MODES,''''];
INITIAL.init_background = BCKGD0;
INITIAL.init_noiselvl = NOISE0;
INITIAL.iseed = 42;
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index badb02a2..eddf3778 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -108,7 +108,6 @@ fprintf(fid,'/\n');
fprintf(fid,'&INITIAL_CON\n');
fprintf(fid,[' INIT_OPT = ', INITIAL.INIT_OPT,'\n']);
-fprintf(fid,[' ACT_ON_MODES = ', INITIAL.ACT_ON_MODES,'\n']);
fprintf(fid,[' init_background = ', num2str(INITIAL.init_background),'\n']);
fprintf(fid,[' init_noiselvl = ', num2str(INITIAL.init_noiselvl),'\n']);
fprintf(fid,[' iseed = ', num2str(INITIAL.iseed),'\n']);
diff --git a/wk/benchmark_and_scan_scripts/old/CBC_kT_nu_scan.m b/wk/CBC_kT_nu_scan.m
similarity index 64%
rename from wk/benchmark_and_scan_scripts/old/CBC_kT_nu_scan.m
rename to wk/CBC_kT_nu_scan.m
index fca6c554..d851b9c9 100644
--- a/wk/benchmark_and_scan_scripts/old/CBC_kT_nu_scan.m
+++ b/wk/CBC_kT_nu_scan.m
@@ -5,24 +5,24 @@ addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
% EXECNAME = 'gyacomo_debug';
-EXECNAME = 'gyacomo';
+EXECNAME = 'gyacomo23_sp';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%
-SIMID = 'p2_linear'; % Name of the simulation
-RERUN = 0; % rerun if the data does not exist
+SIMID = 'p2_linear_new'; % Name of the simulation
+RERUN = 1; % rerun if the data does not exist
RUN = 1;
-KT_a = [3:0.5:6.5 6.96];
-NU_a = [0 0.01:0.01:0.1 0.2 0.5 1.0];
+KT_a = [3.5 4.0 4.5 5.0 5.5 6.0 6.5 6.96];
+NU_a = [0 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2 0.5];
% KT_a = 3.5;
% NU_a = 0.5;
-P = 16;
-J = 8;
+P = 8;
+J = 4;
% collision setting
-CO = 'DG';
-GKCO = 0; % gyrokinetic operator
+CO = 'LD';
+GKCO = 1; % gyrokinetic operator
COLL_KCUT = 1.75;
% model
-KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+NA = 1; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 1e-4; % electron plasma beta
% background gradients setting
K_N = 2.22; % Density '''
@@ -31,7 +31,7 @@ K_N = 2.22; % Density '''
GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
-DT = 2e-3;
+DT = 1e-2;
TMAX = 30;
kymin = 0.3;
NY = 2;
@@ -59,10 +59,8 @@ for KT = KT_a
K_Ne = K_N; % ele Density '''
K_Ni = K_N; % ion Density gradient drive
%% GRID PARAMETERS
- PMAXE = P; % Hermite basis size of electrons
- JMAXE = J; % Laguerre "
- PMAXI = P; % " ions
- JMAXI = J; % "
+ PMAX = P; % Hermite basis size
+ JMAX = J; % Laguerre "
NX = 8; % real space x-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 2*pi/kymin; % Size of the squared frequency domain
@@ -81,12 +79,11 @@ for KT = KT_a
% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
SHIFT_Y = 0.0;
%% TIME PARMETERS
- SPS0D = 1; % Sampling per time unit for 2D arrays
- SPS2D = -1; % Sampling per time unit for 2D arrays
- SPS3D = 1; % Sampling per time unit for 2D arrays
- SPS5D = 1/2; % Sampling per time unit for 5D arrays
- SPSCP = 0; % Sampling per time unit for checkpoints
- JOB2LOAD= -1;
+ DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+ DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+ DTSAVE3D = 2; % Sampling per time unit for 3D arrays
+ DTSAVE5D = 100; % Sampling per time unit for 5D arrays
+ JOB2LOAD = -1; % Start a new simulation serie
%% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
@@ -106,12 +103,18 @@ for KT = KT_a
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
+ ADIAB_E = (NA==1);
HD_CO = 0.0; % Hyper diffusivity cutoff ratio
MU = 0.0; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
+ HYP_V = 'none';
+ HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+ DMAX = -1;
+ NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+ NMAX = 0;
MU_Z = 1.0; %
MU_P = 0.0; %
MU_J = 0.0; %
@@ -126,43 +129,50 @@ for KT = KT_a
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
% check if data exist to run if no data
- data_ = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
- if RUN && (RERUN || isempty(data_.NU_EVOL) || numel(data_.Ts3D)<10)
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
+ data_ = {};
+ try
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ catch
+ data_.outfilenames = [];
+ end
+ if RUN && (RERUN || isempty(data_.outfilenames))
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
+ % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
end
- if ~isempty(data_.NU_EVOL)
- if numel(data_.Ts3D)>10
- % Load results after trying to run
- filename = [SIMID,'/',PARAMS,'/'];
- LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+ if numel(data_.Ts3D)>10
+ % Load results after trying to run
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
- data_ = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
- % linear growth rate (adapted for 2D zpinch and fluxtube)
- options.TRANGE = [0.5 1]*data_.Ts3D(end);
- options.NPLOTS = 0; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
- options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
+ % linear growth rate (adapted for 2D zpinch and fluxtube)
+ options.TRANGE = [0.5 1]*data_.Ts3D(end);
+ options.NPLOTS = 0; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
+ options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
- [~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window
- [~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ...
- field = 0;
- field_t = 0;
- for ik = 2:NY/2+1
- field = squeeze(sum(abs(data_.PHI),3)); % take the sum over z
- field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
- to_measure = log(field_t(it1:it2));
- tw = data_.Ts3D(it1:it2);
- % gr = polyfit(tw,to_measure,1);
- gr = fit(tw,to_measure,'poly1');
- err= confint(gr);
- g_ky(i,j,ik) = gr.p1;
- g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
- end
- [gmax, ikmax] = max(g_ky(i,j,:));
+ [~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window
+ [~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ...
+ field = 0;
+ field_t = 0;
+ for ik = 2:NY/2+1
+ field = squeeze(sum(abs(data_.PHI),3)); % take the sum over z
+ field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
+ to_measure = log(field_t(it1:it2));
+ tw = double(data_.Ts3D(it1:it2));
+% gr = polyfit(tw,to_measure,1);
+ gr = fit(tw,to_measure,'poly1');
+ err= confint(gr);
+ g_ky(i,j,ik) = gr.p1;
+ g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
+ end
+ [gmax, ikmax] = max(g_ky(i,j,:));
- msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.ky(ikmax)); disp(msg);
- end
+ msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
end
i = i + 1;
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index a968d3db..260cd105 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -46,8 +46,8 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
%Paper 2
% folder = '/misc/gene_results/CBC/KT_6.96_64x32x32x24x12_Nexc_5/';
-% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x8x4_Nexc_5_00/';
-folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x16x8_Nexc_5_00/';
+folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x8x4_Nexc_5_00/';
+% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x16x8_Nexc_5_00/';
% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x32x16_Nexc_5_00/';
% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x32x16_Nexc_5_01/';
@@ -66,18 +66,35 @@ folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x16x8_Nexc_5_00/';
% debug ? shearless
% folder = '/misc/gene_results/CBC/shearless_CBC_128x64x24x24x12_00/';
% folder = '/misc/gene_results/CBC/shearless_CBC_128x64x24x24x12_01/';
+if 0
+%% FULL DATA LOAD (LONG)
gene_data = load_gene_data(folder);
+end
gene_data.FIGDIR = folder;
gene_data = invert_kxky_to_kykx_gene_results(gene_data);
-gene_data.grids.Np = gene_data.grids.Nvp-1;
-gene_data.grids.Nj = gene_data.grids.Nmu-1;
+gene_data.grids.Np = gene_data.grids.Nvp;
+gene_data.grids.Nj = gene_data.grids.Nmu;
gene_data.CODENAME = 'GENE';
gene_data.inputs = gene_data.grids;
gene_data.inputs.Na = 1;
gene_data.paramshort = gene_data.params_string;
+if 0
%% Dashboard (Compilation of main plots of the sim)
dashboard(gene_data);
+end
+if 1
+%% ONLY HEAT FLUX
+nrgfile = 'nrg.dat.h5';
+% nrgfile = 'nrg_1.h5';
+T = h5read([folder,nrgfile],'/nrgions/time');
+Qx = h5read([folder,nrgfile],'/nrgions/Q_es');
+[~,it0] = min(abs(0.25*T(end)-T(end));
+Qavg = mean(Qx(it0:end));
+Qstd = std(Qx(it0:end));
+figure
+plot(data.Ts0D,data.HFLUX_X,'DisplayName',folder(32:48))
+end
%% Separated plot routines
if 0
%% Space time diagramm (fig 11 Ivanov 2020)
@@ -88,8 +105,8 @@ options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag
options.INTERP = 1;
options.NCUT = 4; % Number of cuts for averaging and error estimation
options.RESOLUTION = 256;
-fig = plot_radial_transport_and_spacetime(gene_data,options);
-save_figure(gene_data,fig,'.png')
+plot_radial_transport_and_spacetime(gene_data,options);
+% save_figure(gene_data,fig,'.png')
end
if 0
diff --git a/wk/benchmark_and_scan_scripts/GX_2022_fig_1_1.txt b/wk/benchmark_and_scan_scripts/GX_2022_fig_1_1.txt
new file mode 100644
index 00000000..a7cd5b85
--- /dev/null
+++ b/wk/benchmark_and_scan_scripts/GX_2022_fig_1_1.txt
@@ -0,0 +1,11 @@
+0.06956521739130439, 0.013548387096774195
+0.14021739130434785, 0.052258064516129035
+0.2108695652173913, 0.0903225806451613
+0.2804347826086957, 0.12129032258064516
+0.3532608695652174, 0.1393548387096774
+0.4217391304347826, 0.14516129032258066
+0.4945652173913043, 0.13870967741935483
+0.567391304347826, 0.12129032258064516
+0.6380434782608695, 0.09677419354838711
+0.708695652173913, 0.06580645161290323
+0.7782608695652173, 0.03290322580645161
diff --git a/wk/benchmark_and_scan_scripts/GX_2022_fig_2_1.txt b/wk/benchmark_and_scan_scripts/GX_2022_fig_2_1.txt
new file mode 100644
index 00000000..245cf31d
--- /dev/null
+++ b/wk/benchmark_and_scan_scripts/GX_2022_fig_2_1.txt
@@ -0,0 +1,23 @@
+0.07175398633257407, 0.029629629629629617
+0.14236902050113903, 0.09506172839506172
+0.21070615034168572, 0.1617283950617284
+0.28132118451025057, 0.21851851851851853
+0.35649202733485197, 0.2518518518518519
+0.4225512528473805, 0.26296296296296295
+0.4977220956719818, 0.25555555555555554
+0.5637813211845103, 0.2308641975308642
+0.6343963553530751, 0.1962962962962963
+0.7050113895216401, 0.14814814814814817
+0.7779043280182232, 0.09259259259259262
+0.8530751708428246, 0.07654320987654323
+0.9191343963553531, 0.08395061728395065
+0.9920273348519362, 0.09876543209876543
+1.0626423690205011, 0.11234567901234568
+1.1332574031890659, 0.1271604938271605
+1.2038724373576308, 0.14197530864197533
+1.2744874715261958, 0.1617283950617284
+1.347380410022779, 0.1777777777777778
+1.4225512528473803, 0.20246913580246914
+1.4908883826879271, 0.22345679012345682
+1.5546697038724373, 0.24320987654320989
+1.6298405466970387, 0.26790123456790127
diff --git a/wk/benchmark_and_scan_scripts/GX_2022_fig_4_1.txt b/wk/benchmark_and_scan_scripts/GX_2022_fig_4_1.txt
new file mode 100644
index 00000000..6fddbc63
--- /dev/null
+++ b/wk/benchmark_and_scan_scripts/GX_2022_fig_4_1.txt
@@ -0,0 +1,93 @@
+7.9051383399209385, 0.03838771593090229
+28.985507246376812, 0.03838771593090229
+39.52569169960475, 0.7677543186180422
+44.795783926218746, 2.840690978886755
+50.06587615283263, 21.727447216890596
+57.971014492753625, 7.907869481765834
+71.1462450592885, 4.261036468330136
+92.22661396574438, 7.984644913627639
+110.67193675889325, 5.719769673704416
+126.48221343873519, 6.641074856046064
+142.29249011857706, 7.63915547024952
+163.37285902503294, 8.790786948176583
+179.18313570487481, 6.525911708253357
+189.7233201581028, 6.564299424184263
+202.8985507246377, 5.988483685220729
+229.24901185770744, 7.101727447216891
+252.96442687747043, 7.984644913627639
+274.0447957839262, 8.33013435700576
+297.76021080368906, 7.101727447216891
+316.20553359683794, 8.36852207293666
+334.6508563899868, 5.758157389635315
+363.63636363636357, 6.948176583493282
+384.71673254281956, 7.7543186180422286
+411.0671936758893, 6.410748560460654
+437.41765480895907, 7.216890595009598
+450.5928853754942, 6.410748560460654
+479.5783926218708, 8.675623800383878
+498.0237154150197, 6.142034548944338
+511.1989459815547, 7.4856046065259125
+521.7391304347825, 8.90595009596929
+534.9143610013175, 7.869481765834934
+553.3596837944665, 8.982725527831095
+577.0750988142292, 6.602687140115165
+600.790513833992, 7.7543186180422286
+637.68115942029, 6.180422264875244
+650.8563899868248, 7.4856046065259125
+664.0316205533597, 6.833013435700575
+674.5718050065877, 7.984644913627639
+700.9222661396575, 6.564299424184263
+727.2727272727273, 7.984644913627639
+748.353096179183, 7.063339731285989
+764.163372859025, 7.792706333973129
+790.5138339920948, 6.833013435700575
+816.8642951251647, 10.134357005758158
+827.4044795783927, 7.946257197696738
+848.4848484848483, 9.021113243761997
+861.6600790513835, 7.677543186180424
+877.4703557312253, 8.90595009596929
+898.5507246376812, 8.560460652591171
+924.901185770751, 6.218809980806142
+951.2516469038208, 7.7159309021113245
+967.0619235836627, 7.408829174664108
+980.2371541501975, 6.410748560460654
+1006.5876152832673, 8.17658349328215
+1043.478260869565, 7.216890595009598
+1104.0843214756258, 9.328214971209214
+1130.4347826086955, 7.370441458733206
+1148.8801054018445, 8.522072936660269
+1177.8656126482213, 6.679462571976966
+1214.756258234519, 7.44721689059501
+1251.6469038208168, 6.410748560460654
+1270.0922266139657, 7.4856046065259125
+1306.9828722002635, 6.218809980806142
+1341.2384716732543, 8.214971209213052
+1357.0487483530962, 7.63915547024952
+1375.494071146245, 9.289827255278313
+1396.574440052701, 6.90978886756238
+1422.9249011857708, 7.562380038387715
+1441.3702239789195, 6.948176583493282
+1454.5454545454545, 7.7159309021113245
+1470.3557312252965, 6.641074856046064
+1501.97628458498, 7.063339731285989
+1520.421607378129, 6.641074856046064
+1546.772068511199, 8.522072936660269
+1573.1225296442688, 6.871401151631478
+1591.5678524374177, 6.0268714011516344
+1615.2832674571805, 7.408829174664108
+1633.7285902503293, 7.293666026871401
+1665.3491436100132, 8.099808061420346
+1678.5243741765482, 6.602687140115165
+1686.429512516469, 7.332053742802303
+1702.239789196311, 5.758157389635315
+1733.8603425559945, 8.867562380038388
+1757.5757575757575, 6.756238003838774
+1786.561264822134, 7.946257197696738
+1802.3715415019765, 6.410748560460654
+1841.897233201581, 7.024952015355087
+1860.34255599473, 8.445297504798464
+1894.5981554677205, 7.024952015355087
+1913.0434782608695, 7.562380038387715
+1936.758893280632, 6.0268714011516344
+1965.744400527009, 7.4856046065259125
+1989.459815546772, 6.525911708253357
diff --git a/wk/benchmark_and_scan_scripts/GX_2022_fig_4_2.txt b/wk/benchmark_and_scan_scripts/GX_2022_fig_4_2.txt
new file mode 100644
index 00000000..8259b204
--- /dev/null
+++ b/wk/benchmark_and_scan_scripts/GX_2022_fig_4_2.txt
@@ -0,0 +1,5 @@
+1, 0.07058823529411917
+1.498583569405099, 0.14117647058823657
+2, 3.31764705882353
+2.492917847025496, 7.48235294117647
+3, 12.376470588235291
diff --git a/wk/benchmark_and_scan_scripts/old/fort_00.90 b/wk/benchmark_and_scan_scripts/old/fort_00.90
new file mode 100644
index 00000000..b4c87cc3
--- /dev/null
+++ b/wk/benchmark_and_scan_scripts/old/fort_00.90
@@ -0,0 +1,92 @@
+&BASIC
+ nrun = 100000000
+ dt = 0.02
+ tmax = 30
+ maxruntime = 356400
+ job2load = -1
+/
+&GRID
+ pmax = 60
+ jmax = 30
+ Nx = 8
+ Lx = 7.854
+ Ny = 2
+ Ly = 20.944
+ Nz = 24
+ SG = .false.
+ Nexc = 1
+/
+&GEOMETRY
+ geom = 's-alpha'
+ q0 = 1.4
+ shear = 0.8
+ eps = 0.18
+ kappa = 1
+ delta = 0
+ zeta = 0
+ parallel_bc = 'dirichlet'
+ shift_y = 0
+ Npol = 1
+/
+&OUTPUT_PAR
+ dtsave_0d = 1
+ dtsave_1d = -1
+ dtsave_2d = -1
+ dtsave_3d = 2
+ dtsave_5d = 100
+ write_doubleprecision = .false.
+ write_gamma = .true.
+ write_hf = .true.
+ write_phi = .true.
+ write_Na00 = .true.
+ write_Napj = .false.
+ write_dens = .false.
+ write_temp = .true.
+/
+&MODEL_PAR
+ LINEARITY = 'linear'
+ Na = 1
+ mu_x = 0
+ mu_y = 0
+ N_HD = 4
+ mu_z = 1
+ HYP_V = 'hypcoll'
+ mu_p = 0
+ mu_j = 0
+ nu = 0.1
+ k_gB = 1
+ k_cB = 1
+ lambdaD = 0
+ beta = 0.0001
+ ADIAB_E = .true.
+/
+&CLOSURE_PAR
+ hierarchy_closure='truncation'
+ dmax =-1
+ nonlinear_closure='truncation'
+ nmax =0
+/
+&SPECIES
+ name_ = ions
+ tau_ = 1
+ sigma_ = 1
+ q_ = 1
+ K_N_ = 2.22
+ K_T_ = 3.5
+/
+&COLLISION_PAR
+ collision_model = 'DG'
+ GK_CO = .false.
+ INTERSPECIES = .true.
+ mat_file = '/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/oiCa/null'
+ collision_kcut = 1.75
+/
+&INITIAL_CON
+ INIT_OPT = 'phi'
+ init_background = 0
+ init_noiselvl = 1e-05
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
\ No newline at end of file
diff --git a/wk/benchmark_and_scan_scripts/plot_GX_data.m b/wk/benchmark_and_scan_scripts/plot_GX_data.m
new file mode 100644
index 00000000..092a4e72
--- /dev/null
+++ b/wk/benchmark_and_scan_scripts/plot_GX_data.m
@@ -0,0 +1,33 @@
+figure
+% Linear CBC adiab. e, miller (figure 1 top)
+GX_data_ = load('/home/ahoffman/GX_paper_data/fig_1_1.txt');
+subplot(2,2,1)
+scale_x = 1.0;
+scale_y = 2.777778; % gamma_AH/gamma_NM = R/a
+plot(scale_x*GX_data_(:,1),scale_y*GX_data_(:,2),'--ob','Displayname','GX');
+xlabel('$k_y\rho_i$'); ylabel('$\gamma R/v_{ti}$');
+title('(figure 1 top)');
+% Linear CBC kin. e, miller (figure 2 top left)
+GX_data_ = load('/home/ahoffman/GX_paper_data/fig_2_1.txt');
+subplot(2,2,2)
+scale_x = 1.0;
+scale_y = 2.777778; % gamma_AH/gamma_NM = R/a
+plot(scale_x*GX_data_(:,1),scale_y*GX_data_(:,2),'--ob','Displayname','GX');
+xlabel('$k_y\rho_i$'); ylabel('$\gamma R/v_{ti}$');
+title('(figure 2 top left)');
+% Nonlinear heatflux, adiab. e, miller (figure 4 left)
+GX_data_ = load('/home/ahoffman/GX_paper_data/fig_4_1.txt');
+subplot(2,2,3)
+scale_x = 1/2.777778; % t_AH/t_NM = (a/R)
+scale_y = 2.777778^2;
+plot(scale_x*GX_data_(:,1),scale_y*GX_data_(:,2),'-b','Displayname','GX');
+xlabel('$tv_{ti}/R$'); ylabel('$Q_i/G_{GB}^R$');
+title('(figure 4 right)');
+% Dimits shift, adiab. e, miller (figure 4 right)
+GX_data_ = load('/home/ahoffman/GX_paper_data/fig_4_2.txt');
+subplot(2,2,4)
+scale_x = 2.777778; % wT_AH/wT_NM = a/R
+scale_y = 2.777778^2; % Q_AH/Q_NM = (a/R)^2
+plot(scale_x*GX_data_(:,1),scale_y*GX_data_(:,2),'--ob','Displayname','GX');
+xlabel('$R/L_T$'); ylabel('$Q_i/Q_{GB}^R$');
+title('(figure 4 right)');
\ No newline at end of file
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index ec4d127e..91e6f946 100644
--- a/wk/fast_analysis.m
+++ b/wk/fast_analysis.m
@@ -1,86 +1,64 @@
-% Directory of the code "mypathtogyacomo/gyacomo/"
+gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add
+default_plots_options
% Partition of the computer where the data have to be searched
-% PARTITION = '/misc/gyacomo23_outputs/';
-PARTITION = '/home/ahoffman/gyacomo/';
-%% CBC
-% resdir = 'paper_2_GYAC23/CBC/5x3x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/CBC/7x4x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/CBC/9x5x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/CBC/9x5x192x96x32_dp';
-% resdir = 'paper_2_GYAC23/CBC/11x6x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/CBC/11x6x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/CBC/21x11x128x64x24_dp';
-
-%% tests single vs double precision
-% resdir = 'paper_2_GYAC23/precision_study/5x3x128x64x24';
-% resdir = 'paper_2_GYAC23/precision_study/5x3x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/precision_study/5x3x128x64x24_sp';
-% resdir = 'paper_2_GYAC23/precision_study/5x3x128x64x24_sp_clos_1';
-% resdir = 'paper_2_GYAC23/precision_study/3x2x128x64x24_sp_muz_2.0';
-% resdir = 'paper_2_GYAC23/precision_study/test_3x2x128x64x24_sp_muz_2.0';
-% resdir = 'paper_2_GYAC23/precision_study/3x2x128x64x24_sp_clos_1';
-
-%% Marconi results
-% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/5x3x128x64x24_dp_muz_2.0_muxy_0';
-% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/5x3x128x64x24_dp_SG';
-% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/5x3x128x64x24_dp_muz_2.0_full_NL';
-% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/7x4x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/9x5x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/11x6x128x64x24_dp';
-
-% resdir = 'paper_2_GYAC23/collisionless/CBC/5x3x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/collisionless/CBC/7x4x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/collisionless/CBC/9x5x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/collisionless/CBC/11x6x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/collisionless/CBC/9x5x192x96x32_dp';
-
-% resdir = 'paper_2_GYAC23/collisionless/kT_scan_nu_1e-3/5x3x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/collisionless/kT_scan_nu_1e-3/7x4x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/collisionless/kT_scan_nu_1e-3/9x5x128x64x24_dp';
-
-% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/5x3x128x64x24_dp';
-% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/9x5x128x64x24_dp';
-
-%% low precision 3D ITG
-% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/CBC_3x2x64x48x16_CLOS_1';
-% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/kT_0.0';
-% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/kT_3.0';
-% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/kT_3.5';
-% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/kT_4.0';
-% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/kT_4.5';
-% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/kT_5.3';
-% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/CBC';
-
-% resdir = 'results/paper_2_GYAC23/5x2x64x48x16/kT_3.5';
-% resdir = 'results/paper_2_GYAC23/5x2x64x48x16/kT_4.0';
-% resdir = 'results/paper_2_GYAC23/5x2x64x48x16/kT_4.5';
-% resdir = 'results/paper_2_GYAC23/5x2x64x48x16/kT_5.3';
-% resdir = 'results/paper_2_GYAC23/5x2x64x48x16/CBC';
-
-% resdir = 'results/paper_2_GYAC23/9x2x64x48x16/kT_3.5';
-% resdir = 'results/paper_2_GYAC23/9x2x64x48x16/kT_4.0';
-% resdir = 'results/paper_2_GYAC23/9x2x64x48x16/kT_4.5';
-% resdir = 'results/paper_2_GYAC23/9x2x64x48x16/kT_5.3';
-% resdir = 'results/paper_2_GYAC23/9x2x64x48x16/CBC';
-
-% resdir = 'results/paper_2_GYAC23/11x2x64x48x16/kT_3.5';
-% resdir = 'results/paper_2_GYAC23/11x2x64x48x16/kT_4.0';
-% resdir = 'results/paper_2_GYAC23/11x2x64x48x16/kT_4.5';
-% resdir = 'results/paper_2_GYAC23/11x2x64x48x16/kT_5.3';
-% resdir = 'results/paper_2_GYAC23/11x2x64x48x16/CBC';
-
-%% Box size effect on CBC
-% resdir = 'results/paper_2_GYAC23/7x4x128x64x24/CBC_L120';
-resdir = 'results/paper_2_GYAC23/7x4x128x64x24/CBC_L180';
-
-%% testcases
-% resdir = 'testcases/ITG_zpinch';
-% resdir = 'testcases/zpinch_example/results_trunc';
-% resdir = 'testcases/zpinch_example/results_maxd=2';
-% resdir = 'testcases/DLRA_zpinch/base_case';
-% resdir = 'testcases/DLRA_zpinch/nsv_filter_2';
-% resdir = 'testcases/DLRA_zpinch/nsv_filter_6';
-% resdir = 'testcases/cyclone_example';
+PARTITION = '/misc/gyacomo23_outputs/';
+% PARTITION = '/home/ahoffman/gyacomo/';
+
+%% Tests
+% resdir = 'test_stepon_AA/CBC_stepon_AA';
+% resdir = 'test_stepon_AA/CBC_no_stepon_AA'; % No clear conclusions
+%% CBC benchmark
+% resdir = 'paper_2_GYAC23/CBC/3x2x128x64x24';
+% resdir = 'paper_2_GYAC23/CBC/3x2x128x64x24_nu_5e-2';
+% resdir = 'paper_2_GYAC23/CBC/5x3x128x64x24';
+% resdir = 'paper_2_GYAC23/CBC/7x4x128x64x24';
+% resdir = 'paper_2_GYAC23/CBC/21x6x192x96x24';
+
+%% low precision kT scan
+% resdir = 'paper_2_GYAC23/local_runs/1x2x64x48x16/CBC';
+% resdir = 'paper_2_GYAC23/local_runs/1x2x64x48x16/CBC_dp';
+% resdir = 'paper_2_GYAC23/local_runs/1x2x128x64x16/CBC';
+
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/5x2x64x48x16/kT_3.5';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/5x2x64x48x16/kT_4.0';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/5x2x64x48x16/kT_4.5';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/5x2x64x48x16/kT_5.3';
+% resdir = 'paper_2_GYAC23/local_runs/5x2x64x48x16/CBC';
+% resdir = 'paper_2_GYAC23/local_runs/5x2x64x48x16/CBC_noise_init';
+
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/9x2x64x48x16/kT_3.5';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/9x2x64x48x16/kT_4.0';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/9x2x64x48x16/kT_4.5';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/9x2x64x48x16/kT_5.3';
+% resdir = 'paper_2_GYAC23/local_runs/9x2x64x48x16/CBC';
+
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/11x2x64x48x16/kT_4.0';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/11x2x64x48x16/kT_4.5';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/11x2x64x48x16/kT_5.3';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/11x2x64x48x16/CBC';
+% resdir = 'paper_2_GYAC23/local_runs/11x2x64x48x16/CBC';
+
+%% Collision scan
+% resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_CBC/7x2x64x48x16/nu_0.1';
+% resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_CBC/9x2x64x48x16/nu_0.1';
+% resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_CBC/9x2x64x48x16/nu_0.1';
+% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/9x5x128x64x24/nu_0.5';
+% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/5x3x128x64x24/nu_0.5';
+% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x2x128x64x24/nu_0.5';
+% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x2x128x64x24_Lx200/nu_0.5';
+% resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_kT_5.3/9x2x128x64x24/nu_0.01';
+
+%% kT eff study
+% resdir = 'paper_2_GYAC23/kT_eff_study/1x3x128x64x24_kT_3.0/Lx120';
+% resdir = 'paper_2_GYAC23/kT_eff_study/1x3x128x64x24_kT_3.0/Lx240';
+% resdir = 'paper_2_GYAC23/kT_eff_study/3x3x128x64x24_kT_3.4/Lx120';
+% resdir = 'paper_2_GYAC23/kT_eff_study/3x3x128x64x24_kT_3.4/Lx240';
+resdir = 'paper_2_GYAC23/kT_eff_study/5x3x128x64x24_kT_3.5';
+% resdir = 'paper_2_GYAC23/kT_eff_study/7x3x128x64x24_kT_3.6';
%%
J0 = 00; J1 = 10;
@@ -102,15 +80,17 @@ options.NMVA = 1; % Moving average for time traces
options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
options.INTERP = 0;
options.RESOLUTION = 256;
-fig = plot_radial_transport_and_spacetime(data,options);
+plot_radial_transport_and_spacetime(data,options);
end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
+[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
options.INTERP = 1;
options.POLARPLOT = 0;
options.NAME = '\phi';
+% options.NAME = '\phi^{NZ}';
% options.NAME = '\omega_z';
% options.NAME = 'N_i^{00}';
% options.NAME = 's_{Ey}';
@@ -144,8 +124,8 @@ options.NORMALIZE = 0;
options.NAME = 'N_i^{00}';
% options.NAME = '\phi';
options.PLAN = 'xy';
-options.COMP = 'avg';
-options.TIME = [800 900 1000];
+options.COMP = 1;
+options.TIME = [300];
options.RESOLUTION = 256;
fig = photomaton(data,options);
% save_figure(data,fig)
@@ -155,25 +135,25 @@ if 0
profiler(data)
end
-if 0
+if 1
%% Hermite-Laguerre spectrum
-% [data.Napjz, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Napjz');
-[data.Napjz, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'Nipjz');
-options.ST = 0;
+[data.Napjz, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Napjz');
+% [data.Napjz, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'Nipjz');
+options.ST = 1;
options.NORMALIZED = 0;
options.LOGSCALE = 1;
options.FILTER = 0; %filter the 50% time-average of the spectrum from
fig = show_moments_spectrum(data,options);
end
-if 0
+if 1
%% Mode evolution
[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
options.NORMALIZED = 0;
options.TIME = [000:9000];
-options.KX_TW = [1 80]; %kx Growth rate time window
-options.KY_TW = [0 80]; %ky Growth rate time window
+options.KX_TW = [30 40]; %kx Growth rate time window
+options.KY_TW = [10 20]; %ky Growth rate time window
options.NMA = 1;
options.NMODES = 800;
options.iz = 'avg'; % avg or index
diff --git a/wk/fast_heat_flux_analysis.m b/wk/fast_heat_flux_analysis.m
new file mode 100644
index 00000000..6e602d76
--- /dev/null
+++ b/wk/fast_heat_flux_analysis.m
@@ -0,0 +1,43 @@
+%% Analysis of sequential kT scans
+% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/collisionless/kT_scan_nu_1e-3/5x3x128x64x24_dp/';
+% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/collisionless/kT_scan_nu_1e-3/7x4x128x64x24_dp/';
+resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/collisionless/kT_scan_nu_1e-3/9x5x128x64x24_dp/';
+Jobs = [0 2 3 4 5];
+kTA = 0*Jobs;
+kNA = 0*Jobs;
+QxstdA = kTA;
+QxavgA = kTA;
+
+figure
+
+for ij = 1:numel(Jobs)
+ data = {};
+ data = compile_results_low_mem(data,resdir,Jobs(ij),Jobs(ij));
+ % fast heat flux analysis
+ T = data.Ts0D-data.Ts0D(1);
+ Qx = data.HFLUX_X;
+ [~,it0] = min(abs(0.25*T(end)-T));
+ Qavg = mean(Qx(it0:end));
+ Qstd = std(Qx(it0:end));
+ kTA(ij) = data.inputs.K_T;
+ kNA(ij) = data.inputs.K_N;
+ QxavgA(ij) = Qavg;
+ QxstdA(ij) = Qstd;
+ subplot(121)
+ plot(T,Qx,'DisplayName',['$\kappa_T=',num2str(data.inputs.K_T),'$']); hold on
+ plot([T(it0) T(end)],Qavg*[1 1],'--k','DisplayName',...
+ ['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
+end
+[~,indices] = sort(kTA);
+kTA = kTA(indices);
+kNA = kNA(indices);
+QxavgA = QxavgA(indices);
+QxstdA = QxstdA(indices);
+
+% Plots
+
+legend('show')
+subplot(122)
+errorbar(kTA,QxavgA./kTA./kNA,QxstdA./kTA./kNA,'--sr',...
+ 'DisplayName','GYAC LD CBC ($\nu_{DGDK}=0.05$)'); hold on
+xlabel('$\kappa_T$'); ylabel('$\chi$')
diff --git a/wk/heat_flux_convergence_analysis.m b/wk/heat_flux_convergence_analysis.m
new file mode 100644
index 00000000..b80aa6ac
--- /dev/null
+++ b/wk/heat_flux_convergence_analysis.m
@@ -0,0 +1,156 @@
+if 0
+ %% GYAC Local low res
+ % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/nu=0.05/3x2x64x48x16/CBC/';
+ % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/3x2x64x48x16/CBC/';
+ % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/5x2x64x48x16/CBC/';
+ resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/7x3x64x48x16/CBC/';
+ % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/9x2x64x48x16/CBC/';
+ % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/11x2x64x48x16/CBC/';
+ % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/17x2x64x48x16/CBC/';
+ data = {};
+ data = compile_results_low_mem(data,resdir,00,10);
+ % fast heat flux analysis
+ T = data.Ts0D;
+ Qx = data.HFLUX_X;
+ [~,it0] = min(abs(0.25*T(end)-T));
+ Qavg = mean(Qx(it0:end));
+ Qstd = std(Qx(it0:end));
+ figure
+ plot(T,Qx); hold on
+ plot([T(it0) T(end)],Qavg*[1 1],'--k','DisplayName',...
+ ['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
+ legend('show')
+ disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
+end
+if 1
+%% Manually gathered data
+QvsPJ = [...
+ %vp mu Qxav Qxer kT
+ 03 02 62.19 26.46 6.96;...
+ 05 02 51.49 08.38 6.96;...
+ 07 03 50.08 08.38 6.96;...
+ 09 02 42.49 08.52 6.96;...
+ 11 02 36.84 07.22 6.96;... % 192x96 instead of 128x64
+ 17 02 37.26 07.63 6.96;...
+];
+figure
+errorbar(QvsPJ(:,1).*QvsPJ(:,2),QvsPJ(:,3),QvsPJ(:,4),'--s',...
+ 'DisplayName','GYAC 64x48x16 ($\nu_{DGDK}=0.001$)'); hold on
+end
+if 0
+%% Manually gathered data
+QvsPJ = [... % Old results done with nu = 0.05
+ %vp mu Qxav Qxer kT
+ 03 02 97.38 17.42 6.96;...
+ 05 02 64.80 13.50 6.96;...
+ 09 02 54.19 10.56 6.96;...
+ 11 02 59.57 13.58 6.96;... % 192x96 instead of 128x64
+];
+errorbar(QvsPJ(:,1).*QvsPJ(:,2),QvsPJ(:,3),QvsPJ(:,4),'--s',...
+ 'DisplayName','GYAC 64x48x16 ($\nu_{DGDK}=0.05$)'); hold on
+end
+%%%%%%%%%%%%%%%%%%%%%%%%
+if 0
+ %% GYAC Marconi standard res
+% Marconi standard res
+ % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/3x2x128x64x24/';
+% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/3x2x128x64x24_nu_5e-2/';
+% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/5x3x128x64x24/';
+% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/7x4x128x64x24/';
+% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/7x4x128x64x24/CBC_L120/';
+% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/7x4x128x64x24/CBC_L180/';
+% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/21x6x192x96x24/';
+ data = {};
+ data = compile_results_low_mem(data,resdir,00,10);
+ % fast heat flux analysis
+ T = data.Ts0D;
+ Qx = data.HFLUX_X;
+ [~,it0] = min(abs(0.25*T(end)-T));
+ Qavg = mean(Qx(it0:end));
+ Qstd = std(Qx(it0:end));
+ figure
+ plot(T,Qx); hold on
+ plot([T(it0) T(end)],Qavg*[1 1],'--k','DisplayName',...
+ ['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
+ legend('show')
+ disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
+end
+if 1
+%% Manually gathered data
+QvsPJ = [...
+ %vp mu Qxav Qxer kT
+ 03 02 43.60 11.47 6.96;...
+ 03 02 65.03 17.79 6.96;... % nuDGDK = 0.05 instead of 0.001
+ 05 03 44.08 06.51 6.96;...
+ 07 04 33.09 04.85 6.96;...
+ 07 04 35.31 04.15 6.96;... % other run
+ 07 04 37.60 04.65 6.96;... % L=180
+ 21 06 33.07 05.31 6.96;... % 192x96 instead of 128x64
+ % 21 11 00.00 00.00 6.96;... % 192x96 instead of 128x64
+];
+% figure
+errorbar(QvsPJ(:,1).*QvsPJ(:,2),QvsPJ(:,3),QvsPJ(:,4),'--s',...
+ 'DisplayName','GYAC 128x64x24 ($\nu_{DGDK}=0.001$)')
+end
+
+%%%%%%%%%%%%%%%%%%%%%%
+if 0
+ %% GENE
+ % folder = '/misc/gene_results/CBC/KT_6.96_64x32x32x24x12_Nexc_5/';
+% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x8x4_Nexc_5_00/';
+% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x16x8_Nexc_5_00/';
+% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x32x16_Nexc_5_00/';
+% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x32x16_Nexc_5_01/';
+
+ % folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x32x16_Nexc_5_00/';
+ % folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x32x16_Nexc_5_01/';
+ % folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x16x8_Nexc_5/';
+ % folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x8x4_Nexc_5/';
+ folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x8x4_Nexc_5_smallvbox/';
+ % folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x16x8_Nexc_5_largexbox/';
+ % folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x16x8_Muz_0.02/';
+
+
+ % folder = '/misc/gene_results/CBC/128x64x16x6x4/';
+ % fast heat flux analysis
+ nrgfile = 'nrg.dat.h5';
+ % nrgfile = 'nrg_1.h5';
+ T = h5read([folder,nrgfile],'/nrgions/time');
+ Qx = h5read([folder,nrgfile],'/nrgions/Q_es');
+ [~,it0] = min(abs(0.25*T(end)-T));
+ Qavg = mean(Qx(it0:end));
+ Qstd = std(Qx(it0:end));
+ figure
+ plot(T,Qx); hold on
+ plot([T(it0) T(end)],Qavg*[1 1],'--k','DisplayName',...
+ ['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
+ legend('show')
+ disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
+end
+if 1
+%% Manually gathered data
+ QvsNv = [...
+ %vp mu Qxav Qxer kT
+ 06 04 17.90 6.85 6.96;...
+ 08 04 02.01 0.59 6.96;...
+ 08 04 13.12 3.20 6.96;...
+ 16 08 38.13 6.08 6.96;...
+ 16 08 35.74 4.85 6.96;...
+ 32 16 33.10 7.01 6.96;...
+ ];
+% figure
+ errorbar((QvsNv(:,1).*QvsNv(:,2)),QvsNv(:,3),QvsNv(:,4),'--sk',...
+ 'DisplayName','GENE 128x64x24')
+
+% QvsNv = [...
+% %vp mu Qxav Qxer kT
+% 32 16 00.32 0.09 5.30;...
+% 16 08 00.77 0.18 5.30;...
+% ];
+% figure
+% errorbar(QvsNv(:,1).*QvsNv(:,2)/2,QvsNv(:,3),QvsNv(:,4),'sk')
+% legend('GENE KT=5.3');
+end
+set(gca,'xscale','log');
+xlabel('$N_{v\parallel}\times N_{\mu}$'); ylabel('$Q_x$');
+top_title('collisionless CBC');
\ No newline at end of file
diff --git a/wk/kTcrit_convergence_analysis.m b/wk/kTcrit_convergence_analysis.m
new file mode 100644
index 00000000..fe674244
--- /dev/null
+++ b/wk/kTcrit_convergence_analysis.m
@@ -0,0 +1,146 @@
+%% Results of a scan to check the critical kT
+% We set kTcrit as the value where the heat flux can increase by two order
+% of magnitude on 50 t.u. (i.e. from 1e-20 to 1e-18) with an phi backg init
+% of amplitude 1e-5 and on the mode ky=0.2. When the last heat flux value
+% is larger than 1e-15 we fix it as kTmin. Then we decrease kT by 0.1 and
+% rerun. If 1e-15 is not reached, this value is kTcrit. Otherwise iterate.
+% Note: this criterium is equivalent to gamma<0.05
+% P J kTcrit
+kTcritP0 = [...
+0 1 3.5;...
+0 2 2.5;...
+0 3 2.3;...
+0 4 2.4;...
+0 5 2.5;...
+0 6 2.7;...
+0 7 2.9;...
+0 8 2.9;...
+];
+kTcritP1 = [... %used ky=0.3
+1 1 4.4;...
+1 2 3.5;...
+1 3 3.5;...
+1 4 3.8;...
+1 5 4.4;...
+% 1 6 0;...
+% 1 7 0;...
+% 1 8 0;...
+];
+kTcritP2 = [...
+% 2 1 0;...
+% 2 2 0;...
+% 2 3 0;...
+% 2 4 0;...
+% 2 5 0;...
+% 2 6 0;...
+% 2 7 0;...
+% 2 8 0;...
+];
+kTcritP3 = [...
+% 3 1 0;...
+% 3 2 0;...
+% 3 3 0;...
+% 3 4 0;...
+% 3 5 0;...
+% 3 6 0;...
+% 3 7 0;...
+% 3 8 0;...
+];
+kTcritP4 = [...
+4 1 3.5;...
+4 2 2.8;...
+4 3 2.8;...
+4 4 3.0;...
+4 5 4.2;...
+4 6 3.7;...
+4 7 3.4;...
+4 8 3.3;...
+];
+kTcritP5 = [...
+5 1 3.8;...
+5 2 2.9;...
+5 3 3.0;...
+5 4 3.6;...
+5 5 4.2;...
+5 6 3.9;...
+5 7 3.6;...
+5 8 3.5;...
+];
+kTcritP6 = [...
+6 1 3.7;...
+6 2 2.9;...
+6 3 2.9;...
+6 4 3.1;...
+6 5 4.2;...
+6 6 3.8;...
+6 7 3.6;...
+6 8 3.4;...
+];
+kTcritP7 = [...
+7 1 4.0;...
+7 2 3.1;...
+7 3 3.2;...
+7 4 3.7;...
+7 5 4.4;...
+7 6 4.0;...
+7 7 3.7;...
+7 8 3.6;...
+];
+kTcritP8 = [...
+8 1 4.0;...
+8 2 3.1;...
+8 3 3.0;...
+8 4 3.3;...
+8 5 4.4;...
+8 6 4.0;...
+8 7 3.7;...
+8 8 3.5;...
+];
+kTcritP12 = [...
+12 2 3.5;...
+12 4 3.9;...
+12 6 4.3;...
+12 8 3.8;...
+12 10 3.7;...
+12 12 4.0;...
+];
+
+kT_tot = [kTcritP0; kTcritP1; kTcritP2; kTcritP3; kTcritP4; kTcritP5; kTcritP6; kTcritP7; kTcritP8; kTcritP12];
+% Stability map
+Np = 12;
+kT_map = zeros(Np);
+[XX,YY] = meshgrid(0:Np,0:Np);
+sz = size(kT_tot);
+for i = 1:sz(1)
+ kT_map(kT_tot(i,1)+1,kT_tot(i,2)+1) = kT_tot(i,3)-4;
+end
+%plots
+figure
+subplot(121)
+plot(kTcritP0(:,2),kTcritP0(:,3),'--o','DisplayName',num2str(kTcritP0(1,1))); hold on;
+plot(kTcritP4(:,2),kTcritP4(:,3),'--o','DisplayName',num2str(kTcritP4(1,1))); hold on;
+plot(kTcritP5(:,2),kTcritP5(:,3),'--o','DisplayName',num2str(kTcritP5(1,1))); hold on;
+plot(kTcritP6(:,2),kTcritP6(:,3),'--o','DisplayName',num2str(kTcritP6(1,1)));
+plot(kTcritP7(:,2),kTcritP7(:,3),'--o','DisplayName',num2str(kTcritP7(1,1)));
+plot(kTcritP8(:,2),kTcritP8(:,3),'--o','DisplayName',num2str(kTcritP8(1,1)));
+plot(kTcritP12(:,2),kTcritP12(:,3),'--o','DisplayName',num2str(kTcritP12(1,1)));
+ylim([0 5]);
+ylabel('$\kappa_T^{crit}$')
+xlabel('$J$');
+subplot(122)
+imagesc_custom(XX,YY,kT_map);
+colormap(bluewhitered)
+xlabel('$J$');
+ylabel('$P$');
+% bubblechart(kTcritP0(:,1),kTcritP0(:,2),kTcritP0(:,3)); hold on
+% bubblechart(kTcritP4(:,1),kTcritP4(:,2),kTcritP4(:,3)); hold on
+% bubblechart(kTcritP5(:,1),kTcritP5(:,2),kTcritP5(:,3))
+% bubblechart(kTcritP6(:,1),kTcritP6(:,2),kTcritP6(:,3))
+% bubblechart(kTcritP7(:,1),kTcritP7(:,2),kTcritP7(:,3))
+% bubblechart(kTcritP8(:,1),kTcritP8(:,2),kTcritP8(:,3))
+% bubblechart(kTcritP12(:,1),kTcritP12(:,2),kTcritP12(:,3))
+
+
+
+
+
diff --git a/wk/lin_ITG_salpha.m b/wk/lin_ITG_salpha.m
new file mode 100644
index 00000000..9ac342f1
--- /dev/null
+++ b/wk/lin_ITG_salpha.m
@@ -0,0 +1,162 @@
+%% QUICK RUN SCRIPT
+% This script creates a directory in /results and runs a simulation directly
+% from the Matlab framework. It is meant to run only small problems in linear
+% for benchmarking and debugging purposes since it makes Matlab "busy".
+
+%% Set up the paths for the necessary Matlab modules
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
+addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
+addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
+addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
+
+%% Set simulation parameters
+SIMID = 'lin_ITG'; % Name of the simulation
+RUN = 0; % To run or just to load
+default_plots_options
+EXECNAME = 'gyacomo23_sp'; % single precision
+% EXECNAME = 'gyacomo23_dp'; % double precision
+
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NU = 0.001; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_Ne = 0*2.22; % ele Density
+K_Te = 0*6.96; % ele Temperature
+K_Ni = 2.22; % ion Density gradient drive
+K_Ti = 5.3; % ion Temperature
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+NA = 1; % number of kinetic species
+ADIAB_E = (NA==1); % adiabatic electron model
+BETA = 0.0; % electron plasma beta
+%% Set up grid parameters
+P = 60;
+J = 30;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 8; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
+LX = 2*pi/0.05; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.3; % Size of the squared frequency domain in y direction
+NZ = 24; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+
+%% GEOMETRY
+GEOMETRY= 's-alpha';
+% GEOMETRY= 'miller';
+EPS = 0.18; % inverse aspect ratio
+Q0 = 1.4; % safety factor
+SHEAR = 0.8; % magnetic shear
+KAPPA = 1.0; % elongation
+DELTA = 0.0; % triangularity
+ZETA = 0.0; % squareness
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+
+%% TIME PARAMETERS
+TMAX = 50; % Maximal time unit
+DT = 1e-3; % Time step
+DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+DTSAVE3D = 2; % Sampling per time unit for 3D arrays
+DTSAVE5D = 100; % Sampling per time unit for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 0; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
+
+%% OUTPUTS
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.0; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+MU_Z = 2.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 0.0e-5; % Initial noise amplitude
+BCKGD0 = 1.0e-5; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1000; % Cutoff for collision operator
+
+%%-------------------------------------------------------------------------
+%% RUN
+setup
+% system(['rm fort*.90']);
+% Run linear simulation
+if RUN
+ MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
+% RUN =['time mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
+ RUN =['time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
+% RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
+% RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
+ MVOUT='cd ../../../wk;';
+ system([MVIN,RUN,MVOUT]);
+end
+
+%% Analysis
+% load
+filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
+LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
+FIGDIR = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
+% Load outputs from jobnummin up to jobnummax
+J0 = 0; J1 = 0;
+data = {}; % Initialize data as an empty cell array
+% load grids, inputs, and time traces
+data = compile_results_low_mem(data,LOCALDIR,J0,J1);
+
+if 0
+%% Plot heat flux evolution
+figure
+semilogy(data.Ts0D,data.HFLUX_X);
+xlabel('$tc_s/R$'); ylabel('$Q_x$');
+end
+if 1 % Activate or not
+%% plot mode evolution and growth rates
+% Load phi
+[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
+options.NORMALIZED = 0;
+options.TIME = data.Ts3D;
+ % Time window to measure the growth of kx/ky modes
+options.KX_TW = [0.2 1]*data.Ts3D(end);
+options.KY_TW = [0.2 1]*data.Ts3D(end);
+options.NMA = 1; % Set NMA option to 1
+options.NMODES = 999; % Set how much modes we study
+options.iz = 'avg'; % Compressing z
+options.ik = 1; %
+options.fftz.flag = 0; % Set fftz.flag option to 0
+fig = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
+end
+
+
+
diff --git a/wk/lin_Ivanov.m b/wk/lin_Ivanov.m
index ec997fc0..3fd95e20 100644
--- a/wk/lin_Ivanov.m
+++ b/wk/lin_Ivanov.m
@@ -18,10 +18,11 @@ default_plots_options
% EXECNAME = 'gyacomo23_dp'; % double precision
% To compile single precision gyacomo do 'make clean; make sp' in the /gyacomo folder
EXECNAME = 'gyacomo23_sp'; % single precision
+% EXECNAME = 'gyacomo23_debug'; % single precision
%% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
TAU = 1e-2; % e/i temperature ratio
-NU = 0.05/TAU; % Collision frequency
+NU = 0.1/TAU; % Collision frequency
K_Ne = 0*2.22; % ele Density
K_Te = 0*6.96; % ele Temperature
K_Ni = 0*2.22; % ion Density gradient drive
@@ -67,7 +68,7 @@ JOB2LOAD = -1; % Start a new simulation serie
%% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
-CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+CO = 'IV'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
GKCO = 1; % Gyrokinetic operator
ABCO = 1; % INTERSPECIES collisions
INIT_ZF = 0; % Initialize zero-field quantities
diff --git a/wk/load_metadata_scan.m b/wk/load_metadata_scan.m
index df8ed3ec..8dede8e1 100644
--- a/wk/load_metadata_scan.m
+++ b/wk/load_metadata_scan.m
@@ -1,4 +1,10 @@
% Metadata path
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
%% Scans over KT and NU, keeping ky, CO constant (CBC_kT_nu_scan.m)
% datafname = 'p2_linear/8x24_ky_0.3_P_4_J_2_kT_3_6.96_nu_0_1_DGDK.mat';
% datafname = 'p2_linear/8x24_ky_0.3_P_8_J_4_kT_3_6.96_nu_0_1_DGDK.mat';
@@ -13,7 +19,7 @@
% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_SGGK_P_2_20_KT_6.96.mat';
% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_hypcoll_P_2_30_KT_6.96.mat';
% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_hypcoll_P_2_30_KT_5.3.mat';
- datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_dvpar4_P_2_30_KT_6.96.mat';
+ % datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_dvpar4_P_2_30_KT_6.96.mat';
% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_dvpar4_P_2_10_KT_6.96.mat';
% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_dvpar4_P_6_10_KT_6.96.mat';
% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_dvpar4_P_2_10_KT_5.3.mat';
@@ -22,6 +28,15 @@
% datafname = 'p2_linear/8x24_ky_0.3_J_2_15_P_2_30_kT_5.3_DGDK_0.05.mat';
% datafname = 'p2_linear/8x24_ky_0.3_J_2_15_P_2_20_kT_6.96_DGDK_0.02.mat';
% datafname = 'p2_linear/8x24_ky_0.3_J_2_15_P_2_30_kT_5.3_DGDK_0.025.mat';
+%% new
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_4_J_2_kT_3.5_6.96_nu_0_0.5_IVDK.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_4_J_2_kT_3.5_6.96_nu_0_0.5_DGDK.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_4_J_2_kT_3.5_6.96_nu_0_0.5_DGGK.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_4_J_2_kT_3.5_6.96_nu_0_0.5_SGGK.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_4_J_2_kT_3.5_6.96_nu_0_0.5_LDGK.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_6_J_3_kT_3.5_6.96_nu_0_0.5_DGGK.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_6_J_3_kT_3.5_6.96_nu_0_0.5_SGGK.mat';
+datafname = 'p2_linear_new/8x24_ky_0.3_P_6_J_3_kT_3.5_6.96_nu_0_0.5_LDGK.mat';
%% Load data
fname = ['../results/',datafname];
d = load(fname);
diff --git a/wk/local_run.m b/wk/local_run.m
index 36969a57..1e3a3db6 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -15,33 +15,33 @@ addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
SIMID = 'dbg'; % Name of the simulation
RUN = 1; % To run or just to load
default_plots_options
-% EXECNAME = 'gyacomo23_sp'; % single precision
-EXECNAME = 'gyacomo23_dp'; % double precision
+EXECNAME = 'gyacomo23_sp'; % single precision
+% EXECNAME = 'gyacomo23_dp'; % double precision
%% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
-NU = 0.05; % Collision frequency
+NU = 0.001; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_Ne = 0*2.22; % ele Density
K_Te = 0*6.96; % ele Temperature
K_Ni = 1*2.22; % ion Density gradient drive
-K_Ti = 3.6; % ion Temperature
+K_Ti = 5.0; % ion Temperature
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
NA = 1; % number of kinetic species
ADIAB_E = (NA==1); % adiabatic electron model
BETA = 0.0; % electron plasma beta
%% Set up grid parameters
-P = 6;
-J = 2;%P/2;
+P = 1;
+J = 1;%P/2;
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size
-NX = 16; % real space x-gridpoints
-NY = 8; % real space y-gridpoints
+NX = 4; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
-LY = 2*pi/0.2; % Size of the squared frequency domain in y direction
-NZ = 16; % number of perpendicular planes (parallel grid)
+LY = 2*pi/0.3; % Size of the squared frequency domain in y direction
+NZ = 24; % number of perpendicular planes (parallel grid)
SG = 0; % Staggered z grids option
-NEXC = 0; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
%% GEOMETRY
GEOMETRY= 's-alpha';
@@ -57,12 +57,12 @@ SHIFT_Y = 0.0; % Shift in the periodic BC in z
NPOL = 1; % Number of poloidal turns
%% TIME PARAMETERS
-TMAX = 200; % Maximal time unit
-DT = 2e-2; % Time step
+TMAX = 50; % Maximal time unit
+DT = 1e-2; % Time step
DTSAVE0D = 1; % Sampling per time unit for 0D arrays
DTSAVE2D = -1; % Sampling per time unit for 2D arrays
DTSAVE3D = 1; % Sampling per time unit for 3D arrays
-DTSAVE5D = 20; % Sampling per time unit for 5D arrays
+DTSAVE5D = 100; % Sampling per time unit for 5D arrays
JOB2LOAD = -1; % Start a new simulation serie
%% OPTIONS
@@ -99,12 +99,12 @@ MU_X = MU; % Hyperdiffusivity coefficient in x direction
MU_Y = MU; % Hyperdiffusivity coefficient in y direction
N_HD = 4; % Degree of spatial-hyperdiffusivity
MU_Z = 2.0; % Hyperdiffusivity coefficient in z direction
-HYP_V = 'dvpar4'; % Kinetic-hyperdiffusivity model
+HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
LAMBDAD = 0.0; % Lambda Debye
-NOISE0 = 1.0e-5; % Initial noise amplitude
-BCKGD0 = 0.0; % Initial background
+NOISE0 = 0.0e-5; % Initial noise amplitude
+BCKGD0 = 1.0e-5; % Initial background
k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; % Magnetic curvature strength
COLL_KCUT = 1000; % Cutoff for collision operator
@@ -117,7 +117,7 @@ setup
if RUN
MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
% RUN =['time mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
- RUN =['time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0;'];
+ RUN =['time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
% RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
% RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
MVOUT='cd ../../../wk;';
@@ -135,7 +135,7 @@ data = {}; % Initialize data as an empty cell array
% load grids, inputs, and time traces
data = compile_results_low_mem(data,LOCALDIR,J0,J1);
-if 1
+if 0
%% Plot heat flux evolution
figure
semilogy(data.Ts0D,data.HFLUX_X);
diff --git a/wk/multiple_kTscan_p2_analysis.m b/wk/multiple_kTscan_p2_analysis.m
index 4f557fbd..02604c12 100644
--- a/wk/multiple_kTscan_p2_analysis.m
+++ b/wk/multiple_kTscan_p2_analysis.m
@@ -1,44 +1,61 @@
clrs_ = lines(10);
kN=2.22;
-scantype = 'nuscan';
-% scantype = 'kTscan';
+% resdir = '5x3x128x64x24'; clr_ = clrs_(1,:); CBC_res = [44.08 06.51]; kTthresh = 2.8;
+% resdir = '7x4x128x64x24'; clr_ = clrs_(2,:); CBC_res = [44.08 06.51]; kTthresh = 2.9;
+% resdir = '9x2x128x64x24'; clr_ = clrs_(3,:); CBC_res = [37.60 04.65]; kTthresh = 3.3;
+% resdir = '9x2x128x64x24_Lx200'; clr_ = clrs_(3,:); CBC_res = [37.60 04.65]; kTthresh = 3.3;
+% resdir = '9x5x128x64x24'; clr_ = clrs_(3,:); CBC_res = [37.60 04.65]; kTthresh = 3.3;
+% resdir = '9x5x128x64x24_Lx200'; clr_ = clrs_(3,:); CBC_res = [37.60 04.65]; kTthresh = 3.3;
+% resdir = '13x2x128x64x24'; clr_ = clrs_(4,:); CBC_res = [36.84 07.22]; kTthresh = 4.4;
+resdir = '13x5x128x64x24'; clr_ = clrs_(4,:); CBC_res = [37.60 04.65]; kTthresh = 3.9;
-switch scantype
- case 'kTscan'
-% resdir = 'paper_2_GYAC23/collisionless/kT_scan_nu_1e-3/5x3x128x64x24_dp'; clr_ = clrs_(1,:); xname = '$\kappa_T (\kappa_N=2.22)$';
-% resdir = 'paper_2_GYAC23/collisionless/kT_scan_nu_1e-3/7x4x128x64x24_dp'; clr_ = clrs_(2,:); xname = '$\kappa_T (\kappa_N=2.22)$';
- resdir = 'paper_2_GYAC23/collisionless/kT_scan_nu_1e-3/9x5x128x64x24_dp'; clr_ = clrs_(3,:); xname = '$\kappa_T (\kappa_N=2.22)$';
- case 'nuscan'
-% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/5x3x128x64x24_dp'; clr_ = clrs_(1,:); xname = '$\nu (\kappa_T=5.3,\kappa_N=2.22)$';
-% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/9x5x128x64x24_dp'; clr_ = clrs_(3,:); xname = '$\nu (\kappa_T=5.3,\kappa_N=2.22)$';
-
-% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/5x3x128x64x24_dp'; clr_ = clrs_(1,:); xname = '$\nu (\kappa_T=5.3,\kappa_N=2.22)$';
- resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x128x64x24_dp'; clr_ = clrs_(3,:); xname = '$\nu (\kappa_T=5.3,\kappa_N=2.22)$';
+if 1
+ scantype = 'kTscan';
+ xname = '$\kappa_T (\kappa_N=2.22)$';
+ scanvarname = 'kT';
+ GRAD = '';
+ prefix = ['/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/',resdir,'/'];
+ scanval = [6.96 6.5 6.0 5.5 5.0 4.5 4.0];
+ naming = @(s) sprintf('%1.1f',s);
+ titlename = [GRAD,' $\nu=0$'];
+else
+ CO = 'DGGK'; mrkstyl='v'; clr_ = clrs_(2,:);
+ % CO = 'SGGK'; mrkstyl='s'; clr_ = clrs_(1,:);
+ % CO = 'LDGK'; mrkstyl='d'; clr_ = clrs_(5,:);
+ % GRAD = 'CBC';
+ GRAD = 'kT_5.3';
+ % GRAD = 'kT_4.5';
+ scantype = 'nuscan';
+ xname = ['$\nu_{',CO,'}$ '];
+ titlename = [CO,', ',resdir,', ',GRAD];
+ scanvarname = 'nu';
+ prefix = ['/misc/gyacomo23_outputs/paper_2_GYAC23/collision_study/nu',CO,'_scan_',GRAD,'/',resdir,'/'];
+ scanval = [0.005 0.01 0.02 0.05 0.1 0.2 0.5];
+ naming = @(s) num2str(s);
end
-Njobs = 4;
-
-x = 0*(1:Njobs);
-Qx_avg = 0*(1:Njobs);
-Qx_std = 0*(1:Njobs);
-Chi_avg = 0*(1:Njobs);
-Chi_std = 0*(1:Njobs);
+N = numel(scanval);
+x = 1:N;
+Qx_avg = 1:N;
+Qx_std = 1:N;
+Chi_avg = 1:N;
+Chi_std = 1:N;
+data = {};
figure
-datadir = ['/misc/gyacomo23_outputs/',resdir,'/'];
-for i = 1:Njobs+1
- J0 = i-1; J1 = i-1;
+
+for i = 1:N
+ datadir = [prefix,scanvarname,'_',naming(scanval(i)),'/'];
Nseg = 5;
- data = compile_results_low_mem(data,datadir,J0,J1);
+ data = compile_results_low_mem(data,datadir,00,10);
Trange = data.Ts0D(end)*[0.3 1.0];
%
[~,it0] = min(abs(Trange(1) -data.Ts0D));
[~,it1] = min(abs(Trange(end)-data.Ts0D));
%
if 0
- Qx = data.HFLUX_X(it0:it1);
Qxa_ = 0*(1:Nseg);
for n = 1:Nseg
ntseg = floor((it1-it0)/n);
@@ -61,21 +78,63 @@ for i = 1:Njobs+1
case 'kTscan'
x(i) = data.inputs.K_T;
end
+ subplot(N,2,2*i-1)
+ Qx = data.HFLUX_X;
+ T = data.Ts0D;
+ plot(T,Qx,'DisplayName',[scanvarname,'=',num2str(x(i))],...
+ 'Color',clr_); hold on
+ plot([T(it0) T(end)],Qx_avg(i)*[1 1],'--k','DisplayName',...
+ ['$Q_{avg}=',sprintf('%2.2f',Qx_avg(i)),'\pm',sprintf('%2.2f',Qx_std(i)),'$']);
+ ylabel(['$Q_x, (\kappa_T=',num2str(x(i)),')$']);
+ % legend('show')
+end
+% Add colless CBC results
+switch scantype
+ case 'kTscan'
+ Chi_avg = [CBC_res(1)/2.22/6.96 Chi_avg];
+ Chi_std = [CBC_res(2)/2.22/6.96 Chi_std];
+ x = [6.96 x];
+ case 'nuscan'
+ Chi_avg = [CBC_res(1)/2.22/6.96 Chi_avg];
+ Chi_std = [CBC_res(2)/2.22/6.96 Chi_std];
+ x = [0 x];
end
% plot;
-errorbar(x,Chi_avg,Chi_std,'DisplayName',data.paramshort,'color',clr_); hold on;
+subplot(122)
+errorbar(x,Chi_avg,Chi_std,'DisplayName',data.paramshort,'color',clr_,'Marker',mrkstyl); hold on;
switch scantype
case 'nuscan'
Lin1999 = load('/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/Lin_1999_fig2.txt');
plot(Lin1999(:,1),Lin1999(:,2),'--ok','DisplayName','Lin1999');
+ set(gca,'XScale','log')
+ xlim([min(x)*0.8 max(x)*1.2])
case 'kTscan'
Dim2000 = load('/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/Dimits_2000_fig3_full_no_GF.txt');
plot(Dim2000(:,1),Dim2000(:,2),'ok','DisplayName','Dimits 2000');
+ xline(kTthresh,'DisplayName','$\kappa_T^{crit}$','color',clr_)
+ xline(4.0,'DisplayName','Dimits $\kappa_T^{crit}$','color',[0 0 0])
+ xlim([kTthresh*0.8 max(x)*1.2])
+ %-------------- GENE ---------------
+ kT_Qi_GENE = ...
+ [...
+ 13. 2.7e+2 2.2e+1;...%128x64x16x24x12 kymin=0.02 (large box)
+ 11. 1.9e+2 1.7e+1;...%128x64x16x24x12 kymin=0.02 (large box)
+ 9.0 1.1e+2 4.2e+1;...%128x64x16x24x12 kymin=0.05
+ 7.0 3.5e+1 4.6e+0;...%128x64x16x24x12 kymin=0.05
+ 5.3 9.7e+0 6.8e+0;...%128x64x16x24x12 kymin=0.05
+ 4.5 2.3e-1 5.0e-2;...%128x64x16x24x12 kymin=0.05
+ ];
+ y_ = kT_Qi_GENE (:,2)./kT_Qi_GENE (:,1)./kN;
+ e_ = kT_Qi_GENE (:,3)./kT_Qi_GENE (:,1)./kN;
+ errorbar(kT_Qi_GENE (:,1),y_, e_,...
+ '+-.k','DisplayName','GENE 128x64x16x24x12',...
+ 'MarkerSize',msz,'LineWidth',lwt); hold on
end
% plot(ITG_threshold*[1 1],[0 20],'-.','DisplayName','$\kappa_T^{crit}$',...
% 'color',clr_);
ylabel('$\chi$');
xlabel(xname);
-ylim([0,5]);
+title(titlename)
+% ylim(ylimits);
legend('show');
diff --git a/wk/old scripts/header_dp.m b/wk/old scripts/header_dp.m
new file mode 100644
index 00000000..fd4fde5d
--- /dev/null
+++ b/wk/old scripts/header_dp.m
@@ -0,0 +1,67 @@
+%% CBC
+% resdir = 'paper_2_GYAC23/CBC/5x3x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/CBC/7x4x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/CBC/9x5x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/CBC/9x5x192x96x32_dp';
+% resdir = 'paper_2_GYAC23/CBC/11x6x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/CBC/11x6x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/CBC/21x11x128x64x24_dp';
+
+%% tests single vs double precision
+% resdir = 'paper_2_GYAC23/precision_study/5x3x128x64x24';
+% resdir = 'paper_2_GYAC23/precision_study/5x3x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/precision_study/5x3x128x64x24_sp';
+% resdir = 'paper_2_GYAC23/precision_study/5x3x128x64x24_sp_clos_1';
+% resdir = 'paper_2_GYAC23/precision_study/3x2x128x64x24_sp_muz_2.0';
+% resdir = 'paper_2_GYAC23/precision_study/test_3x2x128x64x24_sp_muz_2.0';
+% resdir = 'paper_2_GYAC23/precision_study/3x2x128x64x24_sp_clos_1';
+
+%% Marconi results
+% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/5x3x128x64x24_dp_muz_2.0_muxy_0';
+% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/5x3x128x64x24_dp_SG';
+% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/5x3x128x64x24_dp_muz_2.0_full_NL';
+% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/7x4x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/9x5x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/11x6x128x64x24_dp';
+
+% resdir = 'paper_2_GYAC23/collisionless/CBC/5x3x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collisionless/CBC/7x4x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collisionless/CBC/9x5x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collisionless/CBC/11x6x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collisionless/CBC/9x5x192x96x32_dp';
+
+% resdir = 'paper_2_GYAC23/collisionless/kT_scan_nu_1e-3/5x3x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collisionless/kT_scan_nu_1e-3/7x4x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collisionless/kT_scan_nu_1e-3/9x5x128x64x24_dp';
+
+% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/5x3x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/9x5x128x64x24_dp';
+
+
+%% low precision 3D ITG
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/3x2x64x48x16/CBC_3x2x64x48x16_CLOS_1';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/3x2x64x48x16/kT_0.0';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/3x2x64x48x16/kT_3.0';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/3x2x64x48x16/kT_3.5';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/3x2x64x48x16/kT_4.0';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/3x2x64x48x16/kT_4.5';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/3x2x64x48x16/kT_5.3';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/3x2x64x48x16/CBC';
+
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/5x2x64x48x16/kT_3.5';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/5x2x64x48x16/kT_4.0';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/5x2x64x48x16/kT_4.5';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/5x2x64x48x16/kT_5.3';
+% resdir = 'paper_2_GYAC23/local_runs/5x2x64x48x16/CBC';
+
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/9x2x64x48x16/kT_3.5';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/9x2x64x48x16/kT_4.0';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/9x2x64x48x16/kT_4.5';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/9x2x64x48x16/kT_5.3';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/9x2x64x48x16/CBC';
+
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/11x2x64x48x16/kT_4.0';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/11x2x64x48x16/kT_4.5';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/11x2x64x48x16/kT_5.3';
+% resdir = 'paper_2_GYAC23/local_runs/nu=0.05/11x2x64x48x16/CBC';
+% resdir = 'paper_2_GYAC23/local_runs/11x2x64x48x16/CBC';
diff --git a/wk/scan_lin_ITG.m b/wk/scan_lin_ITG.m
new file mode 100644
index 00000000..d1a81f90
--- /dev/null
+++ b/wk/scan_lin_ITG.m
@@ -0,0 +1,163 @@
+%% QUICK RUN SCRIPT
+% This script creates a directory in /results and runs a simulation directly
+% from the Matlab framework. It is meant to run only small problems in linear
+% for benchmarking and debugging purposes since it makes Matlab "busy".
+
+%% Set up the paths for the necessary Matlab modules
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
+addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
+addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
+addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
+
+%% Set simulation parameters
+SIMID = 'lin_ITG'; % Name of the simulation
+RUN = 0; % To run or just to load
+default_plots_options
+EXECNAME = 'gyacomo23_sp'; % single precision
+% EXECNAME = 'gyacomo23_dp'; % double precision
+
+for i = 1:numel(KT_A)
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NU = 0.001; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_Ne = 0*2.22; % ele Density
+K_Te = 0*6.96; % ele Temperature
+K_Ni = 2.22; % ion Density gradient drive
+K_Ti = KT_A(i); % ion Temperature
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+NA = 1; % number of kinetic species
+ADIAB_E = (NA==1); % adiabatic electron model
+BETA = 0.0; % electron plasma beta
+%% Set up grid parameters
+P = 60;
+J = 30;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 8; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
+LX = 2*pi/0.05; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.3; % Size of the squared frequency domain in y direction
+NZ = 24; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+
+%% GEOMETRY
+GEOMETRY= 's-alpha';
+% GEOMETRY= 'miller';
+EPS = 0.18; % inverse aspect ratio
+Q0 = 1.4; % safety factor
+SHEAR = 0.8; % magnetic shear
+KAPPA = 1.0; % elongation
+DELTA = 0.0; % triangularity
+ZETA = 0.0; % squareness
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+
+%% TIME PARAMETERS
+TMAX = 50; % Maximal time unit
+DT = 1e-3; % Time step
+DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+DTSAVE3D = 2; % Sampling per time unit for 3D arrays
+DTSAVE5D = 100; % Sampling per time unit for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 0; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
+
+%% OUTPUTS
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.0; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+MU_Z = 2.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 0.0e-5; % Initial noise amplitude
+BCKGD0 = 1.0e-5; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1000; % Cutoff for collision operator
+
+%%-------------------------------------------------------------------------
+%% RUN
+setup
+% system(['rm fort*.90']);
+% Run linear simulation
+if RUN
+ MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
+% RUN =['time mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
+ RUN =['time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
+% RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
+% RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
+ MVOUT='cd ../../../wk;';
+ system([MVIN,RUN,MVOUT]);
+end
+
+%% Analysis
+% load
+filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
+LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
+FIGDIR = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
+% Load outputs from jobnummin up to jobnummax
+J0 = 0; J1 = 0;
+data = {}; % Initialize data as an empty cell array
+% load grids, inputs, and time traces
+data = compile_results_low_mem(data,LOCALDIR,J0,J1);
+
+if 0
+%% Plot heat flux evolution
+figure
+semilogy(data.Ts0D,data.HFLUX_X);
+xlabel('$tc_s/R$'); ylabel('$Q_x$');
+end
+if 1 % Activate or not
+%% plot mode evolution and growth rates
+% Load phi
+[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
+options.NORMALIZED = 0;
+options.TIME = data.Ts3D;
+ % Time window to measure the growth of kx/ky modes
+options.KX_TW = [0.2 1]*data.Ts3D(end);
+options.KY_TW = [0.2 1]*data.Ts3D(end);
+options.NMA = 1; % Set NMA option to 1
+options.NMODES = 999; % Set how much modes we study
+options.iz = 'avg'; % Compressing z
+options.ik = 1; %
+options.fftz.flag = 0; % Set fftz.flag option to 0
+fig = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
+end
+end
+
+
--
GitLab