From 83b63e0906121d1b82a8828243c250c3e1867b7e Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Fri, 5 May 2023 09:25:33 +0200
Subject: [PATCH] save scripts
---
matlab/assemble_cosolver_matrices.m | 4 +-
matlab/load/quick_gene_load.m | 3 +-
matlab/plot/plot_cosol_mat.m | 6 +
wk/CBC_kT_PJ_scan.m | 223 ++++++++++++++++++
wk/Dimits_fig3.m | 10 +-
wk/analysis_gene.m | 17 +-
.../Dimits_2000_fig1_gamma.txt | 44 ++++
wk/fast_analysis.m | 13 +-
wk/heat_flux_convergence_analysis.m | 3 +
wk/lin_ITG_salpha.m | 14 +-
wk/lin_Ivanov.m | 8 +-
wk/load_metadata_scan.m | 33 ++-
wk/local_run.m | 21 +-
wk/multiple_kTscan_p2_analysis.m | 140 -----------
wk/nonlin_kT_scan_analysis.m | 200 ++++++++++++++++
wk/nonlin_nu_scan_analysis.m | 142 +++++++++++
16 files changed, 697 insertions(+), 184 deletions(-)
create mode 100644 wk/CBC_kT_PJ_scan.m
create mode 100644 wk/benchmark_and_scan_scripts/Dimits_2000_fig1_gamma.txt
delete mode 100644 wk/multiple_kTscan_p2_analysis.m
create mode 100644 wk/nonlin_kT_scan_analysis.m
create mode 100644 wk/nonlin_nu_scan_analysis.m
diff --git a/matlab/assemble_cosolver_matrices.m b/matlab/assemble_cosolver_matrices.m
index bf71b20b..6bb81548 100644
--- a/matlab/assemble_cosolver_matrices.m
+++ b/matlab/assemble_cosolver_matrices.m
@@ -1,6 +1,6 @@
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/'];
+matname= 'gk_sugama_tau1e-2_P4_J2_dk_5e-2_km_5.0_NFLR_5';
+matdir = [codir,'SGGK_tau1e-2_P4_J2_dk_0.05_kpm_5.0/matrices/'];
kperp = load([matdir,'kperp.in']);
Nkp = numel(kperp);
outdir = '/home/ahoffman/gyacomo/iCa/';
diff --git a/matlab/load/quick_gene_load.m b/matlab/load/quick_gene_load.m
index ef4d4c2e..fb1077c0 100644
--- a/matlab/load/quick_gene_load.m
+++ b/matlab/load/quick_gene_load.m
@@ -79,7 +79,8 @@ path = '/home/ahoffman/gene/linear_CBC_results/';
% 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_60_nw_30_adiabe.txt';
+fname = 'kT_5.3_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';
diff --git a/matlab/plot/plot_cosol_mat.m b/matlab/plot/plot_cosol_mat.m
index 65e04b31..86364624 100644
--- a/matlab/plot/plot_cosol_mat.m
+++ b/matlab/plot/plot_cosol_mat.m
@@ -110,6 +110,7 @@ mfns = {...
'/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_sugama_tau1e-2_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';...
};
@@ -126,6 +127,7 @@ CONAME_A = {...
'LD 11 7 NFLR 16'; ...
'SG 11 7 NFLR 12, tau 1e-3'; ...
'SG 4 2 NFLR 5, tau 1e-3'; ...
+ 'SG 4 2 NFLR 5, tau 1e-2'; ...
% 'Hacked SG A';...
% 'Hacked SG B';...
};
@@ -136,6 +138,7 @@ TAU_A = [1;...
1;...
1e-3;...
1e-3;...
+ 1e-2;...
];
figure
for j_ = 1:numel(mfns)
@@ -177,6 +180,7 @@ mfns = {...
'/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_sugama_tau1e-2_P4_J2_dk_5e-2_km_5.0_NFLR_5.h5';...
};
CONAME_A = {'SG 20 10';...
'PA 20 10';...
@@ -186,6 +190,7 @@ CONAME_A = {'SG 20 10';...
'LD 11 7 NFLR 16';...
'SG 11 7 NFLR 12, tau 1e-3'; ...
'SG 4 2 NFLR 5, tau 1e-3'; ...
+ 'SG 4 2 NFLR 5, tau 1e-2'; ...
};
TAU_A = [1;...
1;...
@@ -195,6 +200,7 @@ TAU_A = [1;...
1;...
1e-3;...
1e-3;...
+ 1e-2;...
];
grow = {};
puls = {};
diff --git a/wk/CBC_kT_PJ_scan.m b/wk/CBC_kT_PJ_scan.m
new file mode 100644
index 00000000..0b428d57
--- /dev/null
+++ b/wk/CBC_kT_PJ_scan.m
@@ -0,0 +1,223 @@
+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';
+EXECNAME = 'gyacomo23_sp';
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%
+SIMID = 'p2_linear_new'; % Name of the simulation
+RERUN = 0; % rerun if the data does not exist
+RUN = 1;
+KT_a = [3:0.5:6.5 6.96];
+% KT_a = [5.5 6];
+% P_a = [25];
+% P_a = [3:1:29];
+P_a = 2:2:10;
+J_a = floor(P_a/2);
+% collision setting
+CO = 'LD';
+NU = 0.1;
+GKCO = 1; % gyrokinetic operator
+COLL_KCUT = 1.75;
+% model
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 1e-4; % electron plasma beta
+% background gradients setting
+K_N = 2.22; % Density '''
+% Geometry
+% GEOMETRY= 'miller';
+GEOMETRY= 's-alpha';
+SHEAR = 0.8; % magnetic shear
+% time and numerical grid
+DT0 = 1e-2;
+TMAX = 30;
+kymin = 0.3;
+NY = 2;
+% arrays for the result
+g_ky = zeros(numel(KT_a),numel(P_a),NY/2+1);
+g_avg= g_ky*0;
+g_std= g_ky*0;
+% Naming of the collision operator
+if GKCO
+ CONAME = [CO,'GK'];
+else
+ CONAME = [CO,'DK'];
+end
+
+j = 1;
+for P = P_a
+i = 1;
+for KT = KT_a
+ %% PHYSICAL PARAMETERS
+ TAU = 1.0; % e/i temperature ratio
+ % SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+ SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+ K_Te = KT; % ele Temperature '''
+ K_Ti = KT; % ion Temperature '''
+ K_Ne = K_N; % ele Density '''
+ K_Ni = K_N; % ion Density gradient drive
+ %% GRID PARAMETERS
+ J = floor(P/2);
+ DT = DT0/sqrt(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
+ NZ = 24; % number of perpendicular planes (parallel grid)
+ NPOL = 1;
+ SG = 0; % Staggered z grids option
+ NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+ %% GEOMETRY
+ % GEOMETRY= 's-alpha';
+ EPS = 0.18; % inverse aspect ratio
+ Q0 = 1.4; % safety factor
+ KAPPA = 1.0; % elongation
+ DELTA = 0.0; % triangularity
+ ZETA = 0.0; % squareness
+ PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
+ SHIFT_Y = 0.0;
+ %% TIME PARMETERS
+ 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
+ ABCO = 1; % INTERSPECIES collisions
+ INIT_ZF = 0; ZF_AMP = 0.0;
+ CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
+ NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
+ KERN = 0; % Kernel model (0 : GK)
+ INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
+ NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
+ %% OUTPUTS
+ W_DOUBLE = 0;
+ W_GAMMA = 1; W_HF = 1;
+ W_PHI = 1; W_NA00 = 1;
+ W_DENS = 0; W_TEMP = 1;
+ W_NAPJ = 0; W_SAPJ = 0;
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % 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; %
+ LAMBDAD = 0.0;
+ NOISE0 = 1.0e-5; % Init noise amplitude
+ BCKGD0 = 0.0; % Init background
+ k_gB = 1.0;
+ k_cB = 1.0;
+ %% RUN
+ setup
+ % naming
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ % check if data exist to run if no data
+ 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
+ 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
+ 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');
+
+ % 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 = 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_.grids.ky(ikmax)); disp(msg);
+ end
+ end
+
+ i = i + 1;
+end
+j = j + 1;
+end
+
+if 0
+%% Check time evolution
+figure;
+to_measure = log(field_t);
+plot(data_.Ts3D,to_measure); hold on
+plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
+end
+
+%% take max growth rate among z coordinate
+y_ = g_ky(:,:,2);
+e_ = g_std(:,:,2);
+
+%%
+if(numel(KT_a)>1 && numel(P_a)>1)
+%% Save metadata
+ktmin = num2str(min(KT_a)); ktmax = num2str(max(KT_a));
+ pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
+filename = [num2str(NX),'x',num2str(NZ),'_ky_',num2str(kymin),...
+ '_kT_',ktmin,'_',ktmax,...
+ '_P_',pmin,'_',pmax,'_',CONAME,'_',num2str(NU),'.mat'];
+metadata.name = filename;
+metadata.kymin = kymin;
+metadata.title = ['$\nu_{',CONAME,'}=$',num2str(NU),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)];
+metadata.par = [num2str(NX),'x1x',num2str(NZ)];
+metadata.nscan = 2;
+metadata.s1name = '$\kappa_T$';
+metadata.s1 = KT_a;
+metadata.s2name = '$P$, $J=\lfloor P/2 \rfloor$';
+metadata.s2 = P_a;
+metadata.dname = '$\gamma c_s/R$';
+metadata.data = y_;
+metadata.err = e_;
+metadata.date = date;
+save([SIMDIR,filename],'-struct','metadata');
+disp(['saved in ',SIMDIR,filename]);
+clear metadata tosave
+end
diff --git a/wk/Dimits_fig3.m b/wk/Dimits_fig3.m
index 3571b46c..fa2e31d9 100644
--- a/wk/Dimits_fig3.m
+++ b/wk/Dimits_fig3.m
@@ -40,7 +40,7 @@
4.5 1.2e+0 5.4e-1;...%192x64x24x6x4 kymin=0.05 ! Lx is too small... (weird oscillations)
];
%-------------- GENE ---------------
- kT_Qi_GENE = ...
+ kT_Qi_GENE_SR = ...
[...
13. 2.7e+2 2.2e+1;...%128x64x16x24x12 kymin=0.02 (large box)
% 13. 2.0e+2 6.6e+1;...%128x64x16x24x12 kymin=0.05
@@ -51,18 +51,22 @@
5.3 9.7e+0 6.8e+0;...%128x64x16x24x12 kymin=0.05
4.5 2.3e-1 5.0e-2;...%128x64x16x24x12 kymin=0.05
];
+ kT_Qi_GENE_HR = ...
+ [...
+ 5.3 3.8e-1 1.7e-1;...%128x64x16x32x16 Nexc=5 kymin=0.05
+ ];
%% Heat conductivity Xi [Ln/rhoi^2/cs] computed as Xi = Qi/kT/kN
%init
kT_Xi_GM_32 = kT_Qi_GM_32;
kT_Xi_GM_53 = kT_Qi_GM_53;
kT_Xi_GM_HD = kT_Qi_GM_HD;
- kT_Xi_GENE = kT_Qi_GENE;
+ kT_Xi_GENE = kT_Qi_GENE_SR;
%scale
for i = 2:3
kT_Xi_GM_32 (:,i) = kT_Qi_GM_32 (:,i)./kT_Qi_GM_32 (:,1)./kN;
kT_Xi_GM_53 (:,i) = kT_Qi_GM_53 (:,i)./kT_Qi_GM_53 (:,1)./kN;
kT_Xi_GM_HD (:,i) = kT_Qi_GM_HD (:,i)./kT_Qi_GM_HD (:,1)./kN;
- kT_Xi_GENE (:,i) = kT_Qi_GENE (:,i)./kT_Qi_GENE (:,1)./kN;
+ kT_Xi_GENE (:,i) = kT_Qi_GENE_SR (:,i)./kT_Qi_GENE_SR (:,1)./kN;
end
%% Dimits fig 3 data
KT_DIM = [4.0 4.5 5.0 6.0 7.0 9.0 12. 14. 16. 18.];
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 260cd105..66ce890e 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -46,13 +46,14 @@ 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_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/KT_4.5_192x96x24x30x16_00/';
+folder = '/misc/gene_results/CBC/KT_5.0_192x96x24x30x16_00/';
% 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/';
@@ -69,7 +70,6 @@ folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x8x4_Nexc_5_00/';
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;
@@ -78,6 +78,7 @@ gene_data.CODENAME = 'GENE';
gene_data.inputs = gene_data.grids;
gene_data.inputs.Na = 1;
gene_data.paramshort = gene_data.params_string;
+end
if 0
%% Dashboard (Compilation of main plots of the sim)
dashboard(gene_data);
@@ -89,11 +90,15 @@ 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));
+[~,it0] = min(abs(T-0.25*T(end)));
Qavg = mean(Qx(it0:end));
-Qstd = std(Qx(it0:end));
+Qstd = std(Qx(it0:end))/2;
figure
-plot(data.Ts0D,data.HFLUX_X,'DisplayName',folder(32:48))
+plot(T,Qx,'DisplayName',folder(32:48)); hold on;
+plot([T(it0) T(end)],Qavg*[1 1],'-k');
+plot([T(it0) T(end)],(Qavg+Qstd)*[1 1],'--k');
+plot([T(it0) T(end)],(Qavg-Qstd)*[1 1],'--k');
+disp(['Q_avg=',sprintf('%2.2e',Qavg),'+-',sprintf('%2.2e',Qstd)]);
end
%% Separated plot routines
if 0
diff --git a/wk/benchmark_and_scan_scripts/Dimits_2000_fig1_gamma.txt b/wk/benchmark_and_scan_scripts/Dimits_2000_fig1_gamma.txt
new file mode 100644
index 00000000..8ff083dd
--- /dev/null
+++ b/wk/benchmark_and_scan_scripts/Dimits_2000_fig1_gamma.txt
@@ -0,0 +1,44 @@
+0.049594813614262545, 0.019378427787934194
+0.07196110210696921, 0.021023765996343702
+0.09918962722852515, 0.02431444241316269
+0.07876823338735822, 0.050091407678244965
+0.0972447325769854, 0.0473491773308958
+0.09918962722852515, 0.04351005484460696
+0.10307941653160452, 0.04076782449725777
+0.14294975688816858, 0.06215722120658136
+0.14683954619124795, 0.07148080438756857
+0.1769854132901134, 0.08190127970749543
+0.16823338735818477, 0.09945155393053018
+0.19935170178282008, 0.08683729433272395
+0.20032414910858992, 0.09396709323583181
+0.19935170178282008, 0.096709323583181
+0.2110210696920583, 0.10109689213893969
+0.22074554294975685, 0.12193784277879344
+0.25575364667747164, 0.10329067641681902
+0.2518638573743922, 0.11206581352833639
+0.26353322528363043, 0.13126142595978063
+0.2878444084278768, 0.11755027422303474
+0.2975688816855753, 0.11151736745886655
+0.2965964343598055, 0.12084095063985376
+0.30923824959481355, 0.12029250457038393
+0.3257698541329011, 0.12029250457038393
+0.3325769854132901, 0.11919561243144426
+0.35105348460291735, 0.12193784277879344
+0.3568881685575364, 0.11645338208409507
+0.3539708265802269, 0.10822669104204755
+0.3539708265802269, 0.10438756855575869
+0.38217179902755266, 0.11151736745886655
+0.4074554294975688, 0.09945155393053018
+0.420097244732577, 0.10383912248628886
+0.4444084278768233, 0.10054844606946985
+0.45899513776337114, 0.09232175502742232
+0.45510534846029166, 0.07312614259597806
+0.4910858995137763, 0.07915904936014627
+0.49692058346839535, 0.056672760511883
+0.5017828200972447, 0.06215722120658136
+0.5319286871961102, 0.07038391224862889
+0.5251215559157212, 0.0627056672760512
+0.5504051863857373, 0.044606946983546614
+0.5990275526742301, 0.023765996343692863
+0.5951377633711508, 0.01992687385740402
+0.6038897893030792, 0.006215722120658129
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index 91e6f946..3da0a52b 100644
--- a/wk/fast_analysis.m
+++ b/wk/fast_analysis.m
@@ -48,8 +48,8 @@ PARTITION = '/misc/gyacomo23_outputs/';
% 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/nuSGGK_scan_kT_5.3/9x5x128x64x24/nu_0.5';
+% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x128x64x24_Lx200/nu_0.5';
% resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_kT_5.3/9x2x128x64x24/nu_0.01';
%% kT eff study
@@ -57,8 +57,9 @@ PARTITION = '/misc/gyacomo23_outputs/';
% 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';
+% resdir = 'paper_2_GYAC23/kT_eff_study/5x3x128x64x24_kT_3.5';
+% resdir = 'paper_2_GYAC23/kT_eff_study/7x3x128x64x24_kT_3.6/truncation';
+% resdir = 'paper_2_GYAC23/kT_eff_study/7x3x128x64x24_kT_3.6/dmax';
%%
J0 = 00; J1 = 10;
@@ -124,8 +125,8 @@ options.NORMALIZE = 0;
options.NAME = 'N_i^{00}';
% options.NAME = '\phi';
options.PLAN = 'xy';
-options.COMP = 1;
-options.TIME = [300];
+options.COMP = 'avg';
+options.TIME = [10 50 80];
options.RESOLUTION = 256;
fig = photomaton(data,options);
% save_figure(data,fig)
diff --git a/wk/heat_flux_convergence_analysis.m b/wk/heat_flux_convergence_analysis.m
index b80aa6ac..80532267 100644
--- a/wk/heat_flux_convergence_analysis.m
+++ b/wk/heat_flux_convergence_analysis.m
@@ -31,7 +31,10 @@ QvsPJ = [...
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
+ 13 02 36.74 08.14 6.96;...
17 02 37.26 07.63 6.96;...
+ 11 06 27.47 05.48 6.96;...
+ 13 05 16.45 04.63 6.96;...
];
figure
errorbar(QvsPJ(:,1).*QvsPJ(:,2),QvsPJ(:,3),QvsPJ(:,4),'--s',...
diff --git a/wk/lin_ITG_salpha.m b/wk/lin_ITG_salpha.m
index 9ac342f1..0df72f58 100644
--- a/wk/lin_ITG_salpha.m
+++ b/wk/lin_ITG_salpha.m
@@ -13,7 +13,7 @@ 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
+RUN = 1; % To run or just to load
default_plots_options
EXECNAME = 'gyacomo23_sp'; % single precision
% EXECNAME = 'gyacomo23_dp'; % double precision
@@ -25,20 +25,20 @@ 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
+K_Ti = 6.96; % 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;
+P = 12;
+J = 6;%P/2;
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size
NX = 8; % real space x-gridpoints
-NY = 2; % real space y-gridpoints
+NY = 12; % 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
+LY = 2*pi/0.1; % 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))
@@ -58,7 +58,7 @@ NPOL = 1; % Number of poloidal turns
%% TIME PARAMETERS
TMAX = 50; % Maximal time unit
-DT = 1e-3; % Time step
+DT = 10e-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
diff --git a/wk/lin_Ivanov.m b/wk/lin_Ivanov.m
index 3fd95e20..505fe95e 100644
--- a/wk/lin_Ivanov.m
+++ b/wk/lin_Ivanov.m
@@ -21,7 +21,7 @@ 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
+TAU = 1e-3; % e/i temperature ratio
NU = 0.1/TAU; % Collision frequency
K_Ne = 0*2.22; % ele Density
K_Te = 0*6.96; % ele Temperature
@@ -68,8 +68,8 @@ JOB2LOAD = -1; % Start a new simulation serie
%% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
-CO = 'IV'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-GKCO = 1; % Gyrokinetic operator
+CO = 'SG'; % 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 = 'max_degree'; % Closure model for higher order moments
@@ -109,7 +109,7 @@ NOISE0 = 1.0e-5; % Initial noise amplitude
BCKGD0 = 0.0; % Initial background
k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; % Magnetic curvature strength
-COLL_KCUT = 1000; % Cutoff for collision operator
+COLL_KCUT = 1.0; % Cutoff for collision operator
%%-------------------------------------------------------------------------
%% RUN
diff --git a/wk/load_metadata_scan.m b/wk/load_metadata_scan.m
index 8dede8e1..728438f4 100644
--- a/wk/load_metadata_scan.m
+++ b/wk/load_metadata_scan.m
@@ -15,9 +15,9 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
% datafname = 'p2_linear/8x24_ky_0.3_kT_3_6.96_P_2_40_DGDK_0.025.mat';
%% Scans over NU and PJ, keeping ky and KY constant (CBC_nu_PJ_scan.m)
% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_DGDK_P_2_40_KT_6.96.mat';
-% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_DGDK_P_2_40_KT_5.3.mat';
+% % datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_DGDK_P_2_40_KT_5.3.mat';
% 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_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_10_KT_6.96.mat';
@@ -34,12 +34,35 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
% 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';
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_6_J_3_kT_3.5_6.96_nu_0_0.5_LDGK.mat';
+
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_8_J_4_kT_3.5_6.96_nu_0_0.5_DGGK.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_8_J_4_kT_3.5_6.96_nu_0_0.5_SGGK.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_8_J_4_kT_3.5_6.96_nu_0_0.5_LDGK.mat';
+
+% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_10_DGDK_0.05.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_30_DGGK_0.05.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_20_SGGK_0.05.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_10_LDGK_0.05.mat';
+
+% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_20_DGDK_0.1.mat';
+datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_20_DGGK_0.1.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_20_SGGK_0.1.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_10_LDGK_0.1.mat';
+
+%% Chose if we filter gamma>0.05
+FILTERGAMMA = 1;
+
%% Load data
fname = ['../results/',datafname];
d = load(fname);
+if FILTERGAMMA
+ d.data = d.data.*(d.data>0.025);
+ d.err = d.err.*(d.data>0.025);
+end
if 1
%% Pcolor of the peak
figure;
@@ -50,8 +73,8 @@ title(d.title);
xlabel(d.s1name); ylabel(d.s2name);
set(gca,'XTicklabel',d.s1)
set(gca,'YTicklabel',d.s2)
-% colormap(jet)
-colormap(bluewhitered)
+colormap(jet)
+% colormap(bluewhitered)
clb=colorbar;
clb.Label.String = '$\gamma c_s/R$';
clb.Label.Interpreter = 'latex';
diff --git a/wk/local_run.m b/wk/local_run.m
index 1e3a3db6..5ba805aa 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -20,25 +20,26 @@ EXECNAME = 'gyacomo23_sp'; % single precision
%% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
-NU = 0.001; % Collision frequency
+NU = 0.05; % 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 = 5.0; % ion Temperature
+K_Ti = 6.96; % 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 = 1;
-J = 1;%P/2;
+P = 16;
+J = P/2;
+DT = 1e-2; % Time step
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size
-NX = 4; % real space x-gridpoints
-NY = 2; % real space y-gridpoints
+NX = 8; % real space x-gridpoints
+NY = 12; % real space y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
-LY = 2*pi/0.3; % Size of the squared frequency domain in y direction
+LY = 2*pi/0.1; % 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))
@@ -58,7 +59,6 @@ NPOL = 1; % Number of poloidal turns
%% TIME PARAMETERS
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
@@ -68,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)
-GKCO = 0; % Gyrokinetic operator
+GKCO = 1; % Gyrokinetic operator
ABCO = 1; % INTERSPECIES collisions
INIT_ZF = 0; % Initialize zero-field quantities
HRCY_CLOS = 'truncation'; % Closure model for higher order moments
@@ -78,6 +78,7 @@ 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)
+% NUMERICAL_SCHEME = 'DOPRI5'; % 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
@@ -141,7 +142,7 @@ figure
semilogy(data.Ts0D,data.HFLUX_X);
xlabel('$tc_s/R$'); ylabel('$Q_x$');
end
-if 0 % Activate or not
+if 1 % Activate or not
%% plot mode evolution and growth rates
% Load phi
[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
diff --git a/wk/multiple_kTscan_p2_analysis.m b/wk/multiple_kTscan_p2_analysis.m
deleted file mode 100644
index 02604c12..00000000
--- a/wk/multiple_kTscan_p2_analysis.m
+++ /dev/null
@@ -1,140 +0,0 @@
-clrs_ = lines(10);
-kN=2.22;
-
-% 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;
-
-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
-
-N = numel(scanval);
-x = 1:N;
-Qx_avg = 1:N;
-Qx_std = 1:N;
-Chi_avg = 1:N;
-Chi_std = 1:N;
-data = {};
-figure
-
-
-for i = 1:N
- datadir = [prefix,scanvarname,'_',naming(scanval(i)),'/'];
- Nseg = 5;
-
- 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
- Qxa_ = 0*(1:Nseg);
- for n = 1:Nseg
- ntseg = floor((it1-it0)/n);
- for m = 1:n
- Qxa_(n) = Qxa_(n) + mean(Qx((1:ntseg)+(m-1)*ntseg));
- end
- Qxa_(n) = Qxa_(n)/n;
- end
- Qx_avg(i) = mean(Qxa_);
- Qx_std(i) = std(Qxa_);
- else
- Qx_avg(i) = mean(data.HFLUX_X(it0:it1));
- Qx_std(i) = std(data.HFLUX_X(it0:it1));
- end
- Chi_avg(i) = Qx_avg(i)./data.inputs.K_T/data.inputs.K_N;
- Chi_std(i) = Qx_std(i)./data.inputs.K_T/data.inputs.K_N;
- switch scantype
- case 'nuscan'
- x(i) = data.inputs.NU;
- 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;
-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);
-title(titlename)
-% ylim(ylimits);
-legend('show');
diff --git a/wk/nonlin_kT_scan_analysis.m b/wk/nonlin_kT_scan_analysis.m
new file mode 100644
index 00000000..1090b797
--- /dev/null
+++ b/wk/nonlin_kT_scan_analysis.m
@@ -0,0 +1,200 @@
+kN=2.22;
+figure
+ERRBAR = 0; LOGSCALE = 0;
+nustr = '1e-3';
+rootdir = ['/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_',nustr];GENE = 0;
+% rootdir = '/misc/gene_results/kT_scan_nu0'; GENE = 1;
+msz = 10; lwt = 2.0;
+mrkstyl='v';
+xname = '$\kappa_T (\kappa_N=2.22)$';
+scanvarname = 'kT';
+scanvalues = [6.96,6.5:-0.5:4.0];
+% Get all subdirectories
+system(['ls -d ',rootdir,'/*/ > list.txt']);
+fid = fopen('list.txt');
+tline = fgetl(fid); i_ = 1; Ps=[]; Js =[]; directories={};
+while ischar(tline)
+ directories{i_} = tline;
+ resstr{i_} = tline(numel(rootdir)+2:end-1);
+ tmp = sscanf(resstr{i_},'%dx%dx%dx%dx%d');
+ Ps = [Ps tmp(1)];
+ Js = [Js tmp(2)];
+ tline = fgetl(fid);
+ i_ = i_+1;
+end
+[~,ids] = sort(Ps);
+fclose(fid);
+system('command rm list.txt');5
+
+directories = directories(ids); Ps = Ps(ids);
+if GENE
+ clrs_ = jet(numel(directories));
+else
+ clrs_ = cool(numel(directories));
+end
+M = numel(directories);
+
+% chi_kT_PJ = zeros(numel(scanvalues),M);
+for j = 1:M
+ % Get all subdirectories
+ system(['ls -d ',directories{j},'*/ > list.txt']);
+ fid = fopen('list.txt');
+ tline = fgetl(fid); i_ = 1; kTs=[]; subdirectories={};
+ while ischar(tline)
+ subdirectories{i_} = tline;
+ str = tline(numel(directories{j})+1:end-1);
+ tmp = sscanf(str,'kT_%f');
+ kTs = [kTs tmp(1)];
+ tline = fgetl(fid);
+ i_ = i_+1;
+ end
+ fclose(fid);
+ system('command rm list.txt');
+ [~,ids] = sort(kTs,'descend');
+ subdirectories = subdirectories(ids); kTs = kTs(ids);
+
+ naming = @(s) sprintf('%1.1f',s); clr_ = clrs_(j,:);
+ N = numel(subdirectories);
+ x = 1:N;
+ Qx_avg = 1:N;
+ Qx_std = 1:N;
+ Chi_avg = 1:N;
+ Chi_std = 1:N;
+ data = {};
+ for i = 1:N
+ subdir = subdirectories{i};
+ if ~GENE
+ data = compile_results_low_mem(data,subdir,00,10);
+ else
+ namelist = read_namelist([subdir,'parameters']);
+ data.inputs.PMAX = namelist.box.nv0;
+ data.inputs.JMAX = namelist.box.nw0;
+ data.inputs.K_T = kTs(i);
+ data.inputs.K_N = namelist.species.omn;
+ nrgfile = 'nrg.dat.h5';
+ data.Ts0D = h5read([subdir,nrgfile],'/nrgions/time');
+ data.HFLUX_X = h5read([subdir,nrgfile],'/nrgions/Q_es');
+ end
+ Trange = data.Ts0D(end)*[0.3 1.0];
+ %
+ [~,it0] = min(abs(Trange(1) -data.Ts0D));
+ [~,it1] = min(abs(Trange(end)-data.Ts0D));
+ %
+ if 0
+ Qxa_ = 0*(1:Nseg);
+ for n = 1:Nseg
+ ntseg = floor((it1-it0)/n);
+ for m = 1:n
+ Qxa_(n) = Qxa_(n) + mean(Qx((1:ntseg)+(m-1)*ntseg));
+ end
+ Qxa_(n) = Qxa_(n)/n;
+ end
+ Qx_avg(i) = mean(Qxa_);
+ Qx_std(i) = std(Qxa_);
+ else
+ Qx_avg(i) = mean(data.HFLUX_X(it0:it1));
+ Qx_std(i) = std(data.HFLUX_X(it0:it1));
+ end
+ Chi_avg(i) = Qx_avg(i)./data.inputs.K_T/data.inputs.K_N;
+ Chi_std(i) = Qx_std(i)./data.inputs.K_T/data.inputs.K_N;
+ x(i) = data.inputs.K_T;
+ subplot(N,2,2*i-1)
+ hold on;
+ Qx = data.HFLUX_X;
+ T = data.Ts0D;
+ plot(T,Qx,'DisplayName',...
+ ['$Q_{avg}=',sprintf('%2.2f',Qx_avg(i)),'\pm',sprintf('%2.2f',Qx_std(i)),'$'],...
+ 'Color',clr_); hold on
+ end
+ % plot;
+ subplot(222)
+ hold on;
+ if ERRBAR
+ errorbar(x,Chi_avg,Chi_std,'DisplayName',...
+ ['(',num2str(data.inputs.PMAX),',',num2str(data.inputs.JMAX),')'],...
+ 'color',clr_,'Marker',mrkstyl);
+ else
+ plot(x,Chi_avg,'DisplayName',...
+ ['(',num2str(data.inputs.PMAX),',',num2str(data.inputs.JMAX),')'],...
+ 'color',clr_,'Marker',mrkstyl);
+ end
+ hold on;
+ % chi_kT_PJ(:,j) = Chi_avg;
+end
+% Formatting
+for i = 1:N
+ subplot(N,2,2*i-1)
+ ylabel('$Q_x$');
+ yl = ylim; xl = xlim;
+ title(['$\kappa_T=',num2str(x(i)),'$'],'Position',[xl(2)/2 yl(2)]);
+ if LOGSCALE
+ set(gca,'YScale','log')
+ else
+ set(gca,'YScale','linear');
+ end
+ if i<N
+ xticklabels([]);
+ else
+ xlabel('$t c_s/R$');
+ end
+end
+subplot(222)
+hold on;
+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(4.0,'DisplayName','Dimits $\kappa_T^{crit}$','color',[0 0 0])
+if LOGSCALE
+ set(gca,'YScale','log')
+end
+ %-------------- 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;
+plot(kT_Qi_GENE (:,1),y_,...
+ '+-.k','DisplayName','GENE 24x12',...
+ 'MarkerSize',msz,'LineWidth',lwt); hold on
+kT_Qi_GENE = ...
+ [...
+ 6.5 1.9e+0;...%128x64x16x32x16 kymin=0.05
+ 6.0 1.0e-1;...%128x64x16x32x16 kymin=0.05
+ 5.5 3.9e-2;...%128x64x16x32x16 kymin=0.05
+ 5.3 1.4e-2;...%128x64x16x32x16 kymin=0.05
+ 4.5 1.2e-2;...%128x64x16x32x16 kymin=0.05
+ 4.0 2.9e-6;...%128x64x16x32x16 kymin=0.05
+ ];
+y_ = kT_Qi_GENE (:,2);
+plot(kT_Qi_GENE (:,1),y_,...
+ '*-.k','DisplayName','GENE 32x16',...
+ 'MarkerSize',msz,'LineWidth',lwt); hold on
+ylabel('$\chi$');
+xlabel(xname);
+title(['$\nu_{DGDK}=$',nustr])
+legend('show');
+legend('Location','northwest')
+xlim([3.5 8]);
+ylim([0 3.5]);
+
+%%
+% subplot(224)
+% figure
+% clrs_ = cool(N);
+% % [PJ,KT] = meshgrid(Ps,scanvalues);
+% % surf(KT,PJ,chi_kT_PJ)
+% for i=1:N
+% target = chi_kT_PJ(i,end);
+% % loglog(Ps,abs(chi_kT_PJ(i,:)-target)/target,'o--','color',clrs_(i,:),...
+% % 'DisplayName',['$\kappa_T=$',num2str(x(i))]); hold on
+% semilogy(Ps,abs(chi_kT_PJ(i,:)-target)/target,'o--','color',clrs_(i,:),...
+% 'DisplayName',['$\kappa_T=$',num2str(x(i))]); hold on
+% end
+% % ylabel('error in \%')
+% ylabel('$\chi$')
+% xlabel('P (J=P/2)');
diff --git a/wk/nonlin_nu_scan_analysis.m b/wk/nonlin_nu_scan_analysis.m
new file mode 100644
index 00000000..5ef4a33f
--- /dev/null
+++ b/wk/nonlin_nu_scan_analysis.m
@@ -0,0 +1,142 @@
+kN=2.22;
+figure
+ERRBAR = 0; LOGSCALE = 0;
+% 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 = '9x5x128x64x24'; 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_Lx200'; clr_ = clrs_(3,:); CBC_res = [37.60 04.65]; kTthresh = 3.3;
+% resdir = '11x6x128x64x24'; clr_ = clrs_(5,:); kTthresh = 3.3;
+% resdir = '13x2x128x64x24'; clr_ = clrs_(4,:); kTthresh = 4.4;
+% resdir = '13x5x128x64x24'; clr_ = clrs_(4,:); kTthresh = 3.9;
+% resdir = '17x9x128x64x24'; clr_ = clrs_(7,:); kTthresh = 3.9;
+resstr={};
+msz = 10; lwt = 2.0;
+% CO = 'DGGK'; mrkstyl='d';
+% CO = 'SGGK'; mrkstyl='s';
+CO = 'LDGK'; mrkstyl='o';
+% GRAD = 'kT_7.0';
+GRAD = 'kT_5.3';
+% GRAD = 'kT_4.5';
+xname = ['$\nu_{',CO,'}$ '];
+titlename = [CO,', ',GRAD];
+scanvarname = 'nu';
+rootdir = ['/misc/gyacomo23_outputs/paper_2_GYAC23/collision_study/nu',CO,'_scan_',GRAD];
+scanval = {'0.005' '0.01' '0.02' '0.05' '0.1' '0.2' '0.5'};
+naming = @(s) num2str(s);
+
+ % Get all directories
+system(['ls -d ',rootdir,'/*/ > list.txt']);
+fid = fopen('list.txt');
+tline = fgetl(fid); i_ = 1; Ps=[]; Js =[]; directories={};
+while ischar(tline)
+ directories{i_} = tline;
+ resstr{i_} = tline(numel(rootdir)+2:end-1);
+ tmp = sscanf(resstr{i_},'%dx%dx%dx%dx%d');
+ Ps = [Ps tmp(1)];
+ Js = [Js tmp(2)];
+ tline = fgetl(fid);
+ i_ = i_+1;
+end
+[~,ids] = sort(Ps);
+fclose(fid);
+system('command rm list.txt');
+
+directories = directories(ids); Ps = Ps(ids);
+clrs_ = cool(numel(directories));
+
+for j = 1:numel(directories)
+ % Get all subdirectories
+ system(['ls -d ',directories{j},'*/ > list.txt']);
+ fid = fopen('list.txt');
+ tline = fgetl(fid); i_ = 1; nus=[]; subdirectories={};
+ while ischar(tline)
+ subdirectories{i_} = tline;
+ str = tline(numel(directories{j})+1:end-1);
+ tmp = sscanf(str,'nu_%f');
+ nus = [nus tmp(1)];
+ tline = fgetl(fid);
+ i_ = i_+1;
+ end
+ fclose(fid);
+ system('command rm list.txt');
+ [~,ids] = sort(nus,'descend');
+ subdirectories = subdirectories(ids); nus = nus(ids);
+
+ naming = @(s) sprintf('%1.1f',s); clr_ = clrs_(j,:);
+ titlename = [CO,', ',GRAD,', ',resstr{j}];
+ N = numel(subdirectories);
+ x = 1:N;
+ Qx_avg = 1:N;
+ Qx_std = 1:N;
+ Chi_avg = 1:N;
+ Chi_std = 1:N;
+ data = {};
+ for i = 1:N
+ subdir = subdirectories{i};
+ data = compile_results_low_mem(data,subdir,00,10);
+ Trange = data.Ts0D(end)*[0.3 1.0];
+ %
+ [~,it0] = min(abs(Trange(1) -data.Ts0D));
+ [~,it1] = min(abs(Trange(end)-data.Ts0D));
+ %
+ Qx_avg(i) = mean(data.HFLUX_X(it0:it1));
+ Qx_std(i) = std(data.HFLUX_X(it0:it1));
+ Chi_avg(i) = Qx_avg(i)./data.inputs.K_T/data.inputs.K_N;
+ Chi_std(i) = Qx_std(i)./data.inputs.K_T/data.inputs.K_N;
+ x(i) = data.inputs.NU;
+ subplot(N,2,2*i-1)
+ hold on;
+ Qx = data.HFLUX_X;
+ T = data.Ts0D;
+ % Plot heatflux vs time
+ 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)),'$']);
+ end
+ % plot chi vs nu
+ subplot(122)
+ hold on;
+ if ERRBAR
+ errorbar(x,Chi_avg,Chi_std,'DisplayName',...
+ ['(',num2str(data.inputs.PMAX),',',num2str(data.inputs.JMAX),')'],...
+ 'color',clr_,'Marker',mrkstyl); hold on;
+ else
+ plot(x,Chi_avg,'DisplayName',...
+ ['(',num2str(data.inputs.PMAX),',',num2str(data.inputs.JMAX),')'],...
+ 'color',clr_,'Marker',mrkstyl);
+ end
+end
+
+% Formatting
+for i = 1:N
+ subplot(N,2,2*i-1)
+ ylabel('$Q_x$');
+ yl = ylim; xl = xlim;
+ yl(1) = 0; ylim(yl);
+ title(['$\nu =',num2str(x(i)),'$'],'Position',[xl(2)/2 yl(2)]);
+ if LOGSCALE
+ set(gca,'YScale','log')
+ else
+ set(gca,'YScale','linear');
+ end
+ if i<N
+ xticklabels([]);
+ else
+ xlabel('$t c_s/R$');
+ end
+end
+
+% ------------- LIN kT=5.3 results
+subplot(122)
+hold on;
+Lin1999 = load('/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/Lin_1999_fig2.txt');
+plot(Lin1999(:,1)/sqrt(2),Lin1999(:,2),'--ok','DisplayName','Lin1999');
+set(gca,'XScale','log')
+xlim([min(x)*0.8 max(x)*1.2])
+ylabel('$\chi$');
+xlabel(xname);
+title(titlename)
+legend('show');
--
GitLab