diff --git a/.gitignore b/.gitignore
index cbb062382e2da7a4de5caef39e883b00da2938fe..da9eb16af31b6cb6fd02ce2ea438592d60848e5e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -39,3 +39,5 @@ local/
*/*.sh
Gallery/
.vscode/settings.json
+*figure*
+*results*
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 08b7fc6c4e5529cb46b4203a401ddbc8de1d89c1..6e1afdce7c2ef6762ad096e1cfc94a730dde0083 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -7,7 +7,7 @@ DATA.NU_EVOL = []; % evolution of parameter nu between jobs
DATA.CO_EVOL = []; % evolution of CO
DATA.MUx_EVOL = []; % evolution of parameter mu between jobs
DATA.MUy_EVOL = []; % evolution of parameter mu between jobs
-DATA.K_Ni_EVOL = []; %
+DATA.K_N_EVOL = []; %
DATA.L_EVOL = []; %
DATA.DT_EVOL = []; %
% FIELDS
@@ -62,7 +62,11 @@ while(CONTINUE)
W_DENS = strcmp(h5readatt(filename,'/data/input','write_dens' ),'y');
W_TEMP = strcmp(h5readatt(filename,'/data/input','write_temp' ),'y');
KIN_E = strcmp(h5readatt(filename,'/data/input', 'KIN_E' ),'y');
- BETA = h5readatt(filename,'/data/input','beta');
+ try
+ BETA = h5readatt(filename,'/data/input','beta');
+ catch
+ BETA = 0;
+ end
% Check polynomials degrees
Pe_new= numel(Pe); Je_new= numel(Je);
Pi_new= numel(Pi); Ji_new= numel(Ji);
@@ -206,7 +210,7 @@ while(CONTINUE)
DATA.CO_EVOL = [DATA.CO_EVOL DATA.CO DATA.CO];
DATA.MUx_EVOL = [DATA.MUx_EVOL DATA.MUx DATA.MUx];
DATA.MUy_EVOL = [DATA.MUy_EVOL DATA.MUy DATA.MUy];
- DATA.K_Ni_EVOL = [DATA.K_Ni_EVOL DATA.K_Ni DATA.K_Ni];
+ DATA.K_N_EVOL = [DATA.K_N_EVOL DATA.K_N DATA.K_N];
DATA.L_EVOL = [DATA.L_EVOL DATA.L DATA.L];
DATA.DT_EVOL = [DATA.DT_EVOL DATA.DT_SIM DATA.DT_SIM];
@@ -281,7 +285,7 @@ else
DATA.dir = DIRECTORY;
DATA.localdir = DIRECTORY;
DATA.param_title=['$\nu_{',DATA.CONAME,'}=$', num2str(DATA.NU), ...
- ', $\kappa_{Ni}=$',num2str(DATA.K_Ni),', $\kappa_{Ti}=$',num2str(DATA.K_Ti),...
+ ', $\kappa_{Ni}=$',num2str(DATA.K_N),', $\kappa_{Ti}=$',num2str(DATA.K_T),...
', $L=',num2str(DATA.L),'$, $N=',...
num2str(DATA.Nx),'$, $(P,J)=(',num2str(DATA.PMAXI),',',...
num2str(DATA.JMAXI),')$,',' $\mu_{hd}=$(',num2str(DATA.MUx),...
diff --git a/matlab/load/load_params.m b/matlab/load/load_params.m
index 89e097229bb61478206863cc04e525f6bb93798f..720e243fc9e6921cdf12a30a14c1fa41952d5855 100644
--- a/matlab/load/load_params.m
+++ b/matlab/load/load_params.m
@@ -1,8 +1,21 @@
DATA.CO = h5readatt(filename,'/data/input','CO');
-DATA.K_Ne = h5readatt(filename,'/data/input','K_Ne');
-DATA.K_Ni = h5readatt(filename,'/data/input','K_Ni');
-DATA.K_Te = h5readatt(filename,'/data/input','K_Te');
-DATA.K_Ti = h5readatt(filename,'/data/input','K_Ti');
+try
+ DATA.ETA_N = h5readatt(filename,'/data/input','ETA_N');
+ DATA.ETA_T = h5readatt(filename,'/data/input','ETA_T');
+catch
+ DATA.ETA_N = 1.0;
+ DATA.ETA_T = 1.0;
+end
+try
+ DATA.K_N = h5readatt(filename,'/data/input','K_n');
+catch
+ DATA.K_N = h5readatt(filename,'/data/input','k_N');
+end
+try
+ DATA.K_T = h5readatt(filename,'/data/input','K_T');
+catch
+ DATA.K_T = h5readatt(filename,'/data/input','k_T');
+end
DATA.Q0 = h5readatt(filename,'/data/input','q0');
DATA.SHEAR = h5readatt(filename,'/data/input','shear');
DATA.EPS = h5readatt(filename,'/data/input','eps');
@@ -21,7 +34,11 @@ DATA.DT_SIM = h5readatt(filename,'/data/input','dt');
DATA.MU = h5readatt(filename,'/data/input','mu');
DATA.MUx = h5readatt(filename,'/data/input','mu_x');
DATA.MUy = h5readatt(filename,'/data/input','mu_y');
-DATA.BETA = h5readatt(filename,'/data/input','beta');
+try
+ DATA.BETA = h5readatt(filename,'/data/input','beta');
+catch
+ DATA.BETA = 0;
+end
DATA.W_GAMMA = h5readatt(filename,'/data/input','write_gamma') == 'y';
DATA.W_PHI = h5readatt(filename,'/data/input','write_phi') == 'y';
DATA.W_NA00 = h5readatt(filename,'/data/input','write_Na00') == 'y';
@@ -48,7 +65,7 @@ else
end
degngrad = [degngrad,'_Kni_%1.1f_nu_%0.0e_',...
DATA.CONAME,'_CLOS_',num2str(DATA.CLOS),'_mu_%0.0e'];
-degngrad = sprintf(degngrad,[DATA.K_Ni,DATA.NU,DATA.MU]);
+degngrad = sprintf(degngrad,[DATA.K_N,DATA.NU,DATA.MU]);
% if ~DATA.LINEARITY; degngrad = ['lin_',degngrad]; end
resolution = [num2str(DATA.Nx),'x',num2str(DATA.Ny),'_'];
gridname = ['L_',num2str(DATA.L),'_'];
diff --git a/matlab/plot/plot_cosol_mat.m b/matlab/plot/plot_cosol_mat.m
index 6940e1d1a07eb10a85eb821c2c7906803189268b..d950808d6b5bf4009144ac1f3e243e302db46360 100644
--- a/matlab/plot/plot_cosol_mat.m
+++ b/matlab/plot/plot_cosol_mat.m
@@ -122,20 +122,23 @@ mfns = {'/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
'/home/ahoffman/HeLaZ/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_4.h5',...
'/home/ahoffman/HeLaZ/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12.h5',...
'/home/ahoffman/HeLaZ/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5',...
+ '/home/ahoffman/HeLaZ/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_30.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
- '/home/ahoffman/HeLaZ/iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5',...
- '/home/ahoffman/HeLaZ/iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5',...
+% '/home/ahoffman/HeLaZ/iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5',...
+% '/home/ahoffman/HeLaZ/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', ...
+ 'FC 10 5 NFLR 30', ...
'FC 4 2 NFLR 6',...
'FC 4 2 NFLR 12', ...
- 'Hacked SG A',...
- 'Hacked SG B'};
+% 'Hacked SG A',...
+% 'Hacked SG B',...
+ };
figure
for j_ = 1:numel(mfns)
mat_file_name = mfns{j_}; disp(mat_file_name);
diff --git a/matlab/setup.m b/matlab/setup.m
index 8600aa5c1d0ff30bd806f35f4a73c62b0bb90f93..f5062e5787b1ab8ea2665e5bc5de8fee4aeb80fc 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -43,10 +43,10 @@ MODEL.q_e =-1.0;
MODEL.q_i = 1.0;
if MODEL.q_e == 0; SIMID = [SIMID,'_i']; end;
% gradients L_perp/L_x
-MODEL.K_Ne = K_Ne; % source term kappa for HW
-MODEL.K_Ni = K_Ni; % source term kappa for HW
-MODEL.K_Te = K_Te; % Temperature
-MODEL.K_Ti = K_Te; % Temperature
+MODEL.K_N = K_N;
+MODEL.ETA_N = ETA_N;
+MODEL.K_T = K_T;
+MODEL.ETA_T = ETA_T;
MODEL.GradB = GRADB; % Magnetic gradient
MODEL.CurvB = CURVB; % Magnetic curvature
MODEL.lambdaD = LAMBDAD;
@@ -99,8 +99,8 @@ else
end
% temp. dens. drives
drives_ = [];
-if abs(K_Ni) > 0; drives_ = [drives_,'_kNi_',num2str(K_Ni)]; end;
-if abs(K_Ti) > 0; drives_ = [drives_,'_kTi_',num2str(K_Ti)]; end;
+if abs(K_N) > 0; drives_ = [drives_,'_kN_',num2str(K_N)]; end;
+if abs(K_T) > 0; drives_ = [drives_,'_kT_',num2str(K_T)]; end;
% collision
coll_ = ['_nu_%1.1e_',CONAME];
coll_ = sprintf(coll_,NU);
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 69fe008dc0338f5d324827ab59ed1dbbde32874e..291fd2e6dfc2ca23714423e81236b01bed823c9c 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -69,10 +69,10 @@ fprintf(fid,[' sigma_e = ', num2str(MODEL.sigma_e),'\n']);
fprintf(fid,[' sigma_i = ', num2str(MODEL.sigma_i),'\n']);
fprintf(fid,[' q_e = ', num2str(MODEL.q_e),'\n']);
fprintf(fid,[' q_i = ', num2str(MODEL.q_i),'\n']);
-fprintf(fid,[' K_Ne = ', num2str(MODEL.K_Ne),'\n']);
-fprintf(fid,[' K_Ni = ', num2str(MODEL.K_Ni),'\n']);
-fprintf(fid,[' K_Te = ', num2str(MODEL.K_Te),'\n']);
-fprintf(fid,[' K_Ti = ', num2str(MODEL.K_Ti),'\n']);
+fprintf(fid,[' K_N = ', num2str(MODEL.K_N),'\n']);
+fprintf(fid,[' ETA_N = ', num2str(MODEL.ETA_N),'\n']);
+fprintf(fid,[' K_T = ', num2str(MODEL.K_T),'\n']);
+fprintf(fid,[' ETA_T = ', num2str(MODEL.ETA_T),'\n']);
fprintf(fid,[' GradB = ', num2str(MODEL.GradB),'\n']);
fprintf(fid,[' CurvB = ', num2str(MODEL.CurvB),'\n']);
fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']);
diff --git a/wk/CBC_kT_PJ_scan.m b/wk/CBC_kT_PJ_scan.m
index 11d27a379b34180c28edea1b68c4e73af6856aba..7d6bae2e193a3ab1e14a47737f448d89b9b7ea1b 100644
--- a/wk/CBC_kT_PJ_scan.m
+++ b/wk/CBC_kT_PJ_scan.m
@@ -3,7 +3,8 @@ default_plots_options
HELAZDIR = '/home/ahoffman/HeLaZ/';
EXECNAME = 'helaz3';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-KT_a = [9:2:17];
+% KT_a = [9:2:17];
+KT_a = 7;
g_max= KT_a*0;
g_avg= KT_a*0;
g_std= KT_a*0;
@@ -15,10 +16,10 @@ DT = 1e-2;
TMAX = 25;
ky_ = 0.3;
SIMID = 'linear_CBC_kT_scan_ky_0.3'; % Name of the simulation
-RUN = 0;
+RUN = 1;
figure
-% P = 12;
-for P = [2 4 6]
+P = 4;
+% for P = [2 4 6]
J = P/2;
i=1;
@@ -114,5 +115,5 @@ hold on;
title(['Linear CBC $K_T$ threshold $k_y=$',num2str(ky_),' (CLOS = 1)']);
legend('show'); xlabel('$K_T$'); ylabel('$\gamma$');
drawnow
-end
+
diff --git a/wk/Dimits_2000_fig3.m b/wk/Dimits_2000_fig3.m
deleted file mode 100644
index 13501cd86b7b315a4a4fef5585ac43d6347076c8..0000000000000000000000000000000000000000
--- a/wk/Dimits_2000_fig3.m
+++ /dev/null
@@ -1,113 +0,0 @@
-%% Heat flux Qi [R/rhos^2/cs]
-kN = 2.22;
-%-------------- GM ---------------
-%(P,J)=(2,1)
-kT_Qi_GM_32 = ...
- [...
- 13. 1.5e+2 1.6e+1;...%192x96x16x3x2 kymin=0.02
- 11. 5.1e+2 3.5e+2;...%192x96x16x3x2 kymin=0.02
- 9.0 9.6e+1 3.0e+1;...%192x96x16x3x2 kymin=0.05
- 7.0 5.0e+1 6.6e+0;...%192x96x16x3x2 kymin=0.05
- 6.0 3.0e+1 4.8e+0;...%192x96x16x3x2 kymin=0.05
- 5.0 1.1e+1 9.4e-1;...%192x96x16x3x2 kymin=0.05
- 4.5 9.2e+0 1.6e+0;...%192x96x16x3x2 kymin=0.05
- ];
-%(P,J)=(4,2)
-kT_Qi_GM_53 = ...
- [...
- 13. 2.0e+2 1.2e+1;...%96x64x16x3x2 kymin=0.02 (large box)
-% 13. 1.1e+2 2.0e+1;...%128x64x16x5x3 kymin=0.02 (large box)
-% 13. 1.3e+2 3.5e+1;...%128x64x16x5x3 kymin=0.05
- 11. 1.2e+2 1.6e+1;...%96x64x16x3x2 kymin=0.02 (large box)
-% 11. 1.6e+2 1.8e+1;...%128x64x16x5x3 kymin=0.02
-% 11. 9.7e+1 2.2e+1;...%128x64x16x5x3 kymin=0.05
- 9.0 8.3e+1 2.2e+1;...%128x64x16x5x3 kymin=0.05
-% 9.0 7.6e+1 2.3e+1;...%256x128x16x3x2 kymin=0.05 (high res)
- 7.0 4.6e+1 2.3e+0;...%128x64x16x5x3 kymin=0.05
- 6.0 3.7e+1 6.9e+0;...%128x64x16x5x3 kymin=0.05
- 5.3 1.9e+1 2.0e+0;...%128x64x16x5x3 kymin=0.05
- 5.0 1.3e+1 3.3e+0;...%128x64x16x5x3 kymin=0.05
- 4.5 9.3e+0 1.0e+0;...%128x64x16x5x3 kymin=0.05
- ];
-%(P,J)=(12,2) or higher
-kT_Qi_GM_122 = ...
- [...
- 7.0 4.1e+1 6.6e+0;...%192x96x24x13x7 kymin=0.05
- 4.5 9.6e-1 1.5e-1;...%128x64x16x13x2 kymin=0.05
- ];
-%-------------- GENE ---------------
-kT_Qi_GENE = ...
- [...
- 13. 2.7e+2 2.2e+1;...%128x64x16x24x12 kymin=0.02 (large box)
-% 13. 2.0e+2 6.6e+1;...%128x64x16x24x12 kymin=0.05
- 11. 1.9e+2 1.7e+1;...%128x64x16x24x12 kymin=0.02 (large box)
-% 11. 3.3e+2 1.6e+2;...%128x64x16x24x12 kymin=0.05
- 9.0 1.1e+2 4.2e+1;...%128x64x16x24x12 kymin=0.05
- 7.0 4.1e+1 2.1e+1;...%128x64x16x24x12 kymin=0.05
- 5.3 1.1e+1 1.8e+1;...%128x64x16x24x12 kymin=0.05
- 4.5 1.9e-1 3.0e-2;...%128x64x16x24x12 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_122 = kT_Qi_GM_122;
-kT_Xi_GENE = kT_Qi_GENE;
-%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_122(:,i) = kT_Qi_GM_122(:,i)./kT_Qi_GM_122(:,1)./kN;
-kT_Xi_GENE (:,i) = kT_Qi_GENE (:,i)./kT_Qi_GENE (:,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.];
-LLNL_GK_DIM = [5.0 0.0; ...
- 7.0 2.5;...
- 9.0 5.0;...
- 12. 8.0;...
- 14. 9.0;...
- 16. 9.5;...
- 18. 10.];
-UCOL_GK_DIM = [5.0 0.5;...
- 6.0 1.0;...
- 7.0 1.5];
-GFL__97_DIM = [4.0 4.0;...
- 4.5 5.0;...
- 5.0 6.5;...
- 7.0 8.0;...
- 9.0 11.];
-GFL__98_DIM = [5.0 1.5;...
- 7.0 5.0;...
- 9.0 7.5];
-%% Plot
-msz = 8.0; lwt = 2.0;
-figure;
-if 1
-xye = kT_Xi_GM_32;
-errorbar(xye(:,1), xye(:,2),xye(:,3),'>-','DisplayName','192x96x16x3x2',...
- 'MarkerSize',msz,'LineWidth',lwt); hold on
-xye = kT_Xi_GM_53;
-errorbar(xye(:,1), xye(:,2),xye(:,3),'<-','DisplayName','128x64x16x5x3',...
- 'MarkerSize',msz,'LineWidth',lwt); hold on
-xye = kT_Xi_GM_122;
-errorbar(xye(:,1), xye(:,2),xye(:,3),'^-','DisplayName','128x64x16x13x3',...
- 'MarkerSize',msz,'LineWidth',lwt); hold on
-xye = kT_Xi_GENE;
-errorbar(xye(:,1), xye(:,2),xye(:,3),'+-.k','DisplayName','GENE 128x64x16x24x12',...
- 'MarkerSize',msz,'LineWidth',lwt); hold on
-end
-if 1
- plot(LLNL_GK_DIM(:,1),LLNL_GK_DIM(:,2),'dk--','DisplayName','Dimits GK, LLNL'); hold on
- plot(UCOL_GK_DIM(:,1),UCOL_GK_DIM(:,2),'*k','DisplayName','Dimits PIC, U.COL');
- plot(GFL__97_DIM(:,1),GFL__97_DIM(:,2),'+k','DisplayName','Dimits GFL 97');
- plot(GFL__98_DIM(:,1),GFL__98_DIM(:,2),'sk','DisplayName','Dimits GFL 98');
-
-end
-plot([6.96 6.96],[0 12],'-.r','DisplayName','CBC');
-plot([4.0 4.00],[0 12],'-.b','DisplayName','$\kappa_T^{crit}$');
-xlabel('$\kappa_T$'); ylabel('$\chi_i$[$L_n/\rho_i^2 v_{thi}$]');
-xlim([0 20]);
-ylim([0 12]);
-title('Dimits et al. 2000, Fig. 3');
-legend('show'); grid on;
\ No newline at end of file
diff --git a/wk/ITG_peak_KT_53_k_035_scan.m b/wk/ITG_peak_KT_53_k_035_scan.m
deleted file mode 100644
index 9671fc0fe8d54e9f14220bca4fb29cebd40709a4..0000000000000000000000000000000000000000
--- a/wk/ITG_peak_KT_53_k_035_scan.m
+++ /dev/null
@@ -1,6 +0,0 @@
-P = [ 2 4 8 10 16];
-
-k = 0.35; KN = 2.22; KT = 5.3;
-
-g = [0.22 0.29 0.10 0.04 0.11];
-w = [0.00 0.00 0.00 0.00 0.00];
\ No newline at end of file
diff --git a/wk/KBM.m b/wk/KBM.m
new file mode 100644
index 0000000000000000000000000000000000000000..f5b5215d9c349266b6cb737412e639e3e62a6b49
--- /dev/null
+++ b/wk/KBM.m
@@ -0,0 +1,188 @@
+%% QUICK RUN SCRIPT
+% This script create a directory in /results and run a simulation directly
+% from matlab framework. It is meant to run only small problems in linear
+% for benchmark and debugging purpose since it makes matlab "busy"
+%
+% SIMID = 'test_circular_geom'; % Name of the simulation
+% SIMID = 'linear_CBC'; % Name of the simulation
+SIMID = 'KBM'; % Name of the simulation
+RUN = 1; % To run or just to load
+addpath(genpath('../matlab')) % ... add
+default_plots_options
+HELAZDIR = '/home/ahoffman/HeLaZ/';
+EXECNAME = 'helaz3';
+% EXECNAME = 'helaz3_dbg';
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set Up parameters
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% PHYSICAL PARAMETERS
+NU = 0.05; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_N = 3.0; % ele Density gradient drive
+ETA_N = 3.0/K_N; % ion Density gradient drive
+K_T = 8.0; % Temperature '''
+ETA_T = 4.5/K_T; % Temperature '''
+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)
+KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 0.03; % electron plasma beta
+%% GRID PARAMETERS
+P = 4;
+J = P/2;
+PMAXE = P; % Hermite basis size of electrons
+JMAXE = J; % Laguerre "
+PMAXI = P; % " ions
+JMAXI = J; % "
+NX = 16; % real space x-gridpoints
+NY = 2; % '' y-gridpoints
+LX = 2*pi/0.1; % Size of the squared frequency domain
+LY = 2*pi/0.25; % Size of the squared frequency domain
+NZ = 16; % number of perpendicular planes (parallel grid)
+NPOL = 1;
+SG = 0; % Staggered z grids option
+%% GEOMETRY
+% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
+GEOMETRY= 's-alpha';
+% GEOMETRY= 'circular';
+Q0 = 1.4; % safety factor
+SHEAR = 0.8; % magnetic shear
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+EPS = 0.18; % inverse aspect ratio
+%% TIME PARMETERS
+TMAX = 25; % Maximal time unit
+DT = 5e-3; % Time step
+SPS0D = 1; % Sampling per time unit for 2D arrays
+SPS2D = 0; % Sampling per time unit for 2D arrays
+SPS3D = 1; % Sampling per time unit for 2D arrays
+SPS5D = 1/5; % Sampling per time unit for 5D arrays
+SPSCP = 0; % Sampling per time unit for checkpoints
+JOB2LOAD= -1;
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+% Collision operator
+% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+CO = 'DG';
+GKCO = 0; % gyrokinetic 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
+%% OUTPUTS
+W_DOUBLE = 1;
+W_GAMMA = 1; W_HF = 1;
+W_PHI = 1; W_NA00 = 1;
+W_DENS = 1; W_TEMP = 1;
+W_NAPJ = 1; W_SAPJ = 0;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% unused
+HD_CO = 0.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;
+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
+GRADB = 1.0;
+CURVB = 1.0;
+%%-------------------------------------------------------------------------
+%% RUN
+setup
+% system(['rm fort*.90']);
+% Run linear simulation
+if RUN
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
+end
+
+%% Load results
+%%
+filename = [SIMID,'/',PARAMS,'/'];
+LOCALDIR = [HELAZDIR,'results/',filename,'/'];
+% Load outputs from jobnummin up to jobnummax
+JOBNUMMIN = 00; JOBNUMMAX = 00;
+data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
+
+%% Short analysis
+if 1
+%% linear growth rate (adapted for 2D zpinch and fluxtube)
+trange = [0.5 1]*data.Ts3D(end);
+nplots = 3;
+lg = compute_fluxtube_growth_rate(data,trange,nplots);
+[gmax, kmax] = max(lg.g_ky(:,end));
+[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
+msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
+msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
+end
+
+if 1
+%% Ballooning plot
+options.time_2_plot = [120];
+options.kymodes = [0.1];
+options.normalized = 1;
+% options.field = 'phi';
+fig = plot_ballooning(data,options);
+end
+
+if 0
+%% Hermite-Laguerre spectrum
+% options.TIME = 'avg';
+options.P2J = 1;
+options.ST = 1;
+options.PLOT_TYPE = 'space-time';
+% options.PLOT_TYPE = 'Tavg-1D';
+% options.PLOT_TYPE = 'Tavg-2D';
+options.NORMALIZED = 0;
+options.JOBNUM = 0;
+options.TIME = [0 50];
+options.specie = 'i';
+options.compz = 'avg';
+fig = show_moments_spectrum(data,options);
+% fig = show_napjz(data,options);
+save_figure(data,fig)
+end
+
+if 0
+%% linear growth rate for 3D Zpinch (kz fourier transform)
+trange = [0.5 1]*data.Ts3D(end);
+options.keq0 = 1; % chose to plot planes at k=0 or max
+options.kxky = 1;
+options.kzkx = 0;
+options.kzky = 0;
+[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
+save_figure(data,fig)
+end
+if 0
+%% Mode evolution
+options.NORMALIZED = 0;
+options.K2PLOT = 1;
+options.TIME = [0:1000];
+options.NMA = 1;
+options.NMODES = 1;
+options.iz = 'avg';
+fig = mode_growth_meter(data,options);
+save_figure(data,fig,'.png')
+end
+
+
+if 0
+%% RH TEST
+ikx = 2; t0 = 0; t1 = data.Ts3D(end);
+[~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D));
+plt = @(x) squeeze(mean(real(x(1,ikx,:,it0:it1)),3))./squeeze(mean(real(x(1,ikx,:,it0)),3));
+figure
+plot(data.Ts3D(it0:it1), plt(data.PHI));
+xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
+title(sprintf('$k_x=$%2.2f, $k_y=0.00$',data.kx(ikx)))
+end
diff --git a/wk/RH_test.m b/wk/RH_test.m
new file mode 100644
index 0000000000000000000000000000000000000000..463a73af7eb9d901f65dfa9481655fa90f817e2d
--- /dev/null
+++ b/wk/RH_test.m
@@ -0,0 +1,188 @@
+%% QUICK RUN SCRIPT
+% This script create a directory in /results and run a simulation directly
+% from matlab framework. It is meant to run only small problems in linear
+% for benchmark and debugging purpose since it makes matlab "busy"
+%
+% SIMID = 'test_circular_geom'; % Name of the simulation
+% SIMID = 'linear_CBC'; % Name of the simulation
+SIMID = 'RH_test'; % Name of the simulation
+RUN = 1; % To run or just to load
+addpath(genpath('../matlab')) % ... add
+default_plots_options
+HELAZDIR = '/home/ahoffman/HeLaZ/';
+EXECNAME = 'helaz3';
+% EXECNAME = 'helaz3_dbg';
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set Up parameters
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% PHYSICAL PARAMETERS
+NU = 0.0; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_N = 0.0; % ele Density gradient drive
+ETA_N = 0.0/K_N; % ion Density gradient drive
+K_T = 0.0; % Temperature '''
+ETA_T = 0.0/K_T; % Temperature '''
+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)
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 0.0; % electron plasma beta
+%% GRID PARAMETERS
+P = 20;
+J = P/2;
+PMAXE = P; % Hermite basis size of electrons
+JMAXE = J; % Laguerre "
+PMAXI = P; % " ions
+JMAXI = J; % "
+NX = 2; % real space x-gridpoints
+NY = 2; % '' y-gridpoints
+LX = 2*pi/0.05; % Size of the squared frequency domain
+LY = 2*pi/0.25; % Size of the squared frequency domain
+NZ = 16; % number of perpendicular planes (parallel grid)
+NPOL = 1;
+SG = 0; % Staggered z grids option
+%% GEOMETRY
+% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
+GEOMETRY= 's-alpha';
+% GEOMETRY= 'circular';
+Q0 = 1.4; % safety factor
+SHEAR = 0.0; % magnetic shear
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+EPS = 0.18; % inverse aspect ratio
+%% TIME PARMETERS
+TMAX = 100; % Maximal time unit
+DT = 5e-3; % Time step
+SPS0D = 1; % Sampling per time unit for 2D arrays
+SPS2D = 0; % Sampling per time unit for 2D arrays
+SPS3D = 1; % Sampling per time unit for 2D arrays
+SPS5D = 1/5; % Sampling per time unit for 5D arrays
+SPSCP = 0; % Sampling per time unit for checkpoints
+JOB2LOAD= -1;
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+% Collision operator
+% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+CO = 'DG';
+GKCO = 0; % gyrokinetic 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= 'mom00'; % Start simulation with a noisy mom00/phi/allmom
+%% OUTPUTS
+W_DOUBLE = 1;
+W_GAMMA = 1; W_HF = 1;
+W_PHI = 1; W_NA00 = 1;
+W_DENS = 1; W_TEMP = 1;
+W_NAPJ = 1; W_SAPJ = 0;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% unused
+HD_CO = 0.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;
+MU_Z = 1.0; %
+MU_P = 0.0; %
+MU_J = 0.0; %
+LAMBDAD = 0.0;
+NOISE0 = 1.0e-5; % Init noise amplitude
+BCKGD0 = 1.0; % Init background
+GRADB = 1.0;
+CURVB = 1.0;
+%%-------------------------------------------------------------------------
+%% RUN
+setup
+% system(['rm fort*.90']);
+% Run linear simulation
+if RUN
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
+end
+
+%% Load results
+%%
+filename = [SIMID,'/',PARAMS,'/'];
+LOCALDIR = [HELAZDIR,'results/',filename,'/'];
+% Load outputs from jobnummin up to jobnummax
+JOBNUMMIN = 00; JOBNUMMAX = 00;
+data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
+
+%% Short analysis
+if 0
+%% linear growth rate (adapted for 2D zpinch and fluxtube)
+trange = [0.5 1]*data.Ts3D(end);
+nplots = 3;
+lg = compute_fluxtube_growth_rate(data,trange,nplots);
+[gmax, kmax] = max(lg.g_ky(:,end));
+[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
+msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
+msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
+end
+
+if 0
+%% Ballooning plot
+options.time_2_plot = [120];
+options.kymodes = [0.1];
+options.normalized = 1;
+% options.field = 'phi';
+fig = plot_ballooning(data,options);
+end
+
+if 0
+%% Hermite-Laguerre spectrum
+% options.TIME = 'avg';
+options.P2J = 1;
+options.ST = 1;
+options.PLOT_TYPE = 'space-time';
+% options.PLOT_TYPE = 'Tavg-1D';
+% options.PLOT_TYPE = 'Tavg-2D';
+options.NORMALIZED = 0;
+options.JOBNUM = 0;
+options.TIME = [0 50];
+options.specie = 'i';
+options.compz = 'avg';
+fig = show_moments_spectrum(data,options);
+% fig = show_napjz(data,options);
+save_figure(data,fig)
+end
+
+if 0
+%% linear growth rate for 3D Zpinch (kz fourier transform)
+trange = [0.5 1]*data.Ts3D(end);
+options.keq0 = 1; % chose to plot planes at k=0 or max
+options.kxky = 1;
+options.kzkx = 0;
+options.kzky = 0;
+[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
+save_figure(data,fig)
+end
+if 0
+%% Mode evolution
+options.NORMALIZED = 0;
+options.K2PLOT = 1;
+options.TIME = [0:1000];
+options.NMA = 1;
+options.NMODES = 1;
+options.iz = 'avg';
+fig = mode_growth_meter(data,options);
+save_figure(data,fig,'.png')
+end
+
+
+if 1
+%% RH TEST
+ikx = 2; t0 = 0; t1 = data.Ts3D(end);
+[~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D));
+plt = @(x) squeeze(mean(real(x(1,ikx,:,it0:it1)),3))./squeeze(mean(real(x(1,ikx,:,it0)),3));
+figure
+plot(data.Ts3D(it0:it1), plt(data.PHI));
+xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
+title(sprintf('$k_x=$%2.2f, $k_y=0.00$',data.kx(ikx)))
+end
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index dfa4d3a87c1e4aa5533c8e02f01c94a0edd77d94..fbfbd936d360f9a9174f02915d447e2389cfe3ff 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -9,7 +9,7 @@ MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
system(['mkdir -p ',MISCDIR]);
system(['mkdir -p ',LOCALDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
-system(CMD);
+% system(CMD);
% Load outputs from jobnummin up to jobnummax
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
data.localdir = LOCALDIR;
@@ -23,8 +23,8 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
% data.scale = 1;%/(data.Nx*data.Ny)^2;
-options.TAVG_0 = 25;%0.4*data.Ts3D(end);
-options.TAVG_1 = 40;%0.9*data.Ts3D(end); % Averaging times duration
+options.TAVG_0 = 250;%0.4*data.Ts3D(end);
+options.TAVG_1 = 350;%0.9*data.Ts3D(end); % Averaging times duration
options.NCUT = 4; % Number of cuts for averaging and error estimation
options.NMVA = 1; % Moving average for time traces
% options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
@@ -51,7 +51,7 @@ options.NAME = '\phi';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'kxky';
+options.PLAN = 'xz';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
@@ -105,11 +105,11 @@ if 0
options.SPAR = gene_data.vp';
options.XPERP = gene_data.mu';
options.iz = 'avg';
-options.T = [300 600];
+options.T = [250 600];
options.PLT_FCT = 'contour';
options.ONED = 0;
options.non_adiab = 0;
-options.SPECIE = 'e';
+options.SPECIE = 'i';
options.RMS = 1; % Root mean square i.e. sqrt(sum_k|f_k|^2) as in Gene
fig = plot_fa(data,options);
% save_figure(data,fig,'.png')
@@ -121,14 +121,14 @@ if 0
options.P2J = 0;
options.ST = 1;
options.PLOT_TYPE = 'space-time';
-options.NORMALIZED = 0;
+options.NORMALIZED = 1;
options.JOBNUM = 0;
options.TIME = [1000];
options.specie = 'i';
options.compz = 'avg';
fig = show_moments_spectrum(data,options);
% fig = show_napjz(data,options);
-save_figure(data,fig,'.png');
+% save_figure(data,fig,'.png');
end
if 0
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 62fb9452a8f0e648191e022e0255c4f74dfddab0..0ae250f5c4137f19cdb3013675b175fa99850868 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -15,23 +15,23 @@ addpath(genpath([helazdir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/HP_fig_2a_mu_1e-2/';
% folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
% folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
-folder = '/misc/gene_results/LD_zpinch_1.6/';
+% folder = '/misc/gene_results/LD_zpinch_1.6/';
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
% folder = '/misc/gene_results/Z-pinch/ZP_HP_kn_1.6_HRES/';
% folder = '/misc/gene_results/ZP_kn_2.5_large_box/';
% folder = '/misc/gene_results/CBC/128x64x16x24x12/';
-% folder = '/misc/gene_results/CBC/196x96x20x32x16_01/';
+folder = '/misc/gene_results/CBC/196x96x20x32x16_01/';
% folder = '/misc/gene_results/CBC/128x64x16x6x4/';
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x16x24x12_01/';
-folder = '/misc/gene_results/CBC/KT_4.5_128x64x16x24x12_01/';
-% folder = '/misc/gene_results/CBC/KT_11_128x64x16x24x12/';
-% folder = '/misc/gene_results/CBC/KT_11_large_box_128x64x16x24x12/';
+% folder = '/misc/gene_results/CBC/KT_4.5_128x64x16x24x12_01/';
+% folder = '/misc/gene_results/CBC/KT_9_128x64x16x24x12/';
+% folder = '/misc/gene_results/CBC/KT_13_large_box_128x64x16x24x12/';
gene_data = load_gene_data(folder);
gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-options.TAVG_0 = 0.6*gene_data.Ts3D(end);
+options.TAVG_0 = 0.1*gene_data.Ts3D(end);
options.TAVG_1 = gene_data.Ts3D(end); % Averaging times duration
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, Q_x)
diff --git a/wk/cbc_qx_results.m b/wk/cbc_qx_results.m
deleted file mode 100644
index 2fd7e9b08c20d17a16f4073879a19ff25c1b7a6e..0000000000000000000000000000000000000000
--- a/wk/cbc_qx_results.m
+++ /dev/null
@@ -1,41 +0,0 @@
-%% CBC BENCHMARK
-cbc = [0080 0100 0120];
-gm42 = [15.4 32.2 43.2];
-gm42_err = [2.22 05.2 08.1];
-gm84 = [11.3 25.2 44.7];
-gm84_err = [2.45 03.4 06.2];
-gne = [8.44 24.1 43.5];
-gne_err = [1.55 04.3 05.9];
-kN = [1.776 2.22 2.664];
-kT = [5.568 6.96 8.352];
-figure
-errorbar(cbc,gm42,gm42_err,'o-','LineWidth',1.5); hold on;
-errorbar(cbc,gm84,gm84_err,'o-','LineWidth',1.5);
-errorbar(cbc,gne,gne_err,'x-k','LineWidth',1.5);
-
-% set(gca, 'YScale', 'log')
-
-legend('GM (4,2)','GM (8,4)','Gene')
-xlabel('CBC drive [\%]'); ylabel('Radial Heat Flux $Q_x^\infty$');
-
-
-%% DIMITS
-
-KN = 2.22;
-KT = [1.00 0.90 0.80 0.70 0.60 0.50]*6.96;
-gm42 = [32.2 18.8 10.5 5.89 1.74 0.00];
-gm42_err = [05.2 03.8 2.10 1.66 0.53 0.00];
-gm84 = [0.00 13.2 7.66 2.86 0.00 0.00];
-gm84_err = [0.00 2.79 1.93 0.55 0.00 0.00];
-gne = [0.00 0.00 0.00 0.00 0.00 0.00];
-gne_err = [0.00 0.00 0.00 0.00 0.00 0.00];
-
-figure
-errorbar(KT,gm42./KT,gm42_err./KT,'o-', 'LineWidth',1.5); hold on;
-errorbar(KT,gm84./KT,gm84_err./KT,'o-', 'LineWidth',1.5);
-errorbar(KT, gne./KT, gne_err./KT,'x-k','LineWidth',1.5);
-
-% set(gca, 'YScale', 'log')
-
-legend('GM (4,2)','GM (8,4)','Gene')
-xlabel('$\eta=\kappa_T/\kappa_N$'); ylabel('$\chi_i$');
diff --git a/wk/g2k3_transport_scaling.m b/wk/g2k3_transport_scaling.m
deleted file mode 100644
index c75cd19861a5916a4fe16e244501d07b2be254a3..0000000000000000000000000000000000000000
--- a/wk/g2k3_transport_scaling.m
+++ /dev/null
@@ -1,38 +0,0 @@
-%% G ~ g^2/k^3 rule
-
-K_N = [1.60 1.70 1.80 1.90 2.00 2.10 2.20 2.30 2.40 2.50];
-%DG
-g_DG = [0.08 0.20 0.31 0.42 0.55 0.61 0.70 0.79 0.87 0.96];
-k_DG = [0.21 0.31 0.37 0.42 0.47 0.47 0.47 0.52 0.52 0.47];
-G_DG = g_DG.^2./k_DG.^3./K_N;
-%SG
-g_SG = [0.08 0.19 0.30 0.40 0.49 0.59 0.68 0.76 0.85 0.93];
-k_SG = [0.26 0.37 0.42 0.42 0.47 0.47 0.52 0.52 0.52 0.47];
-G_SG = g_SG.^2./k_SG.^3./K_N;
-%LR
-g_LR = [0.12 0.18 0.25 0.32 0.39 0.47 0.55 0.64 0.72 0.81];
-k_LR = [0.31 0.31 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37];
-G_LR = g_LR.^2./k_LR.^3./K_N;
-%LD
-g_LD = [0.05 0.11 0.19 0.27 0.35 0.44 0.54 0.63 0.73 0.83];
-k_LD = [0.26 0.26 0.31 0.31 0.31 0.31 0.31 0.31 0.31 0.31];
-G_LD = g_LD.^2./k_LD.^3./K_N;
-% nu=0
-g_CL = [0.22 0.27 0.32 0.37 0.44 0.51 0.59 0.67 0.76 0.85];
-k_CL = [0.73 0.73 0.73 0.73 0.63 0.52 0.52 0.52 0.52 0.42];
-G_CL = g_CL.^2./k_CL.^3./K_N;
-%% Plot
-
-figure;
-% plot(K_N,G_SG); hold on;
-% plot(K_N,G_DG); hold on;
-% plot(K_N,G_LR); hold on;
-% plot(K_N,G_LD); hold on;
-% plot(K_N,G_CL,'--k'); hold on;
-k = 0.3;
-plot(K_N,g_SG.^2/k.^3./K_N); hold on;
-plot(K_N,g_DG.^2/k.^3./K_N); hold on;
-plot(K_N,g_LR.^2/k.^3./K_N); hold on;
-plot(K_N,g_LD.^2/k.^3./K_N); hold on;
-plot(K_N,g_CL.^2/k.^3,'--k'); hold on;
-% plot(K_N,(G_DG+G_SG+G_LR+G_LD)/4);
\ No newline at end of file
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 5ccf79448facc6d60186d0f25e437777194955b3..1874ce85fd72ab361cd63bc6f73a9b059939094c 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -50,12 +50,13 @@ helazdir = '/home/ahoffman/HeLaZ/';
% outfile = 'CBC/kT_4.5_192x96x24x13x7';
% outfile = 'CBC/kT_5.3_192x96x24x13x7';
% outfile = 'CBC/kT_13_large_box_128x64x16x5x3';
-outfile = 'CBC/kT_13_96x96x16x3x2_Nexc_6';
+% outfile = 'CBC/kT_13_96x96x16x3x2_Nexc_6';
% outfile = 'CBC/kT_11_96x64x16x5x3_ky_0.02';
% outfile = 'CBC/kT_scan_128x64x16x5x3';
% outfile = 'CBC/kT_scan_192x96x16x3x2';
+outfile = 'CBC/kT_13_96x96x16x3x2_Nexc_6';
%% Linear CBC
% outfile = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe';
diff --git a/wk/quick_run.m b/wk/quick_run.m
index 23acafa52edc071b2d22bef29e55ef0881c648d0..aff73ce3c67ffa7442095134aa1293c202b0dec6 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -6,7 +6,7 @@
% SIMID = 'test_circular_geom'; % Name of the simulation
% SIMID = 'linear_CBC'; % Name of the simulation
SIMID = 'dbg'; % Name of the simulation
-RUN = 1; % To run or just to load
+RUN = 0; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
HELAZDIR = '/home/ahoffman/HeLaZ/';
@@ -19,25 +19,25 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%% PHYSICAL PARAMETERS
NU = 0.05; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_Ne = 3.0; % ele Density gradient drive
-K_Ni = 3.0;%2.0; % ion Density gradient drive
-K_Te = 4.5; % Temperature '''
-K_Ti = 8.0;%0.25*K_N; % Temperature '''
+K_N = 2.22; % ele Density gradient drive
+ETA_N = 1.0; % ion Density gradient drive
+K_T = 6.96; % Temperature '''
+ETA_T = 1.0; % Temperature '''
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)
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 0.0; % electron plasma beta
%% GRID PARAMETERS
-P = 4;
-J = 2;%P/2;
+P = 12;
+J = P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
-NX = 32; % real space x-gridpoints
-NY = 16; % '' y-gridpoints
+NX = 16; % real space x-gridpoints
+NY = 8; % '' y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain
-LY = 2*pi/0.05; % Size of the squared frequency domain
+LY = 2*pi/0.1; % Size of the squared frequency domain
NZ = 16; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
@@ -47,10 +47,10 @@ GEOMETRY= 's-alpha';
% GEOMETRY= 'circular';
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
-NEXC = 6; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 25; % Maximal time unit
+TMAX = 40; % Maximal time unit
DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays