diff --git a/matlab/compute/compute_fluxtube_growth_rate.m b/matlab/compute/compute_fluxtube_growth_rate.m
index 47ba061c751089e261e319abff8578a2a9890ef0..7ee54f667b992b0014f740806c7481ca6c154b36 100644
--- a/matlab/compute/compute_fluxtube_growth_rate.m
+++ b/matlab/compute/compute_fluxtube_growth_rate.m
@@ -6,8 +6,8 @@ if DATA.Nx > 1
else
ikxnz = abs(DATA.PHI(1,:,1,1)) > -1;
end
-ikynz = abs(DATA.PHI(:,1,1,1)) > 0;
-
+ikynz = (abs(DATA.PHI(:,1,1,1)) > 0);
+ikynz = logical(ikynz .* (DATA.ky > 0));
phi = DATA.PHI(ikynz,ikxnz,:,:);
t = DATA.Ts3D;
@@ -39,34 +39,44 @@ linear_gr.g_ky = real(w_ky(Is,:));
linear_gr.w_ky = imag(w_ky(Is,:));
linear_gr.ce = abs(ce);
linear_gr.ky = kys;
+linear_gr.std_g = std (real(w_ky(Is,:)),[],2);
+linear_gr.avg_g = mean(real(w_ky(Is,:)),2);
+linear_gr.std_w = std (imag(w_ky(Is,:)),[],2);
+linear_gr.avg_w = mean(imag(w_ky(Is,:)),2);
+
if PLOT >0
- figure
+ figure
+if PLOT >1
+ subplot(1,PLOT,1)
+end
+
plot(linear_gr.ky,linear_gr.g_ky(:,end),'-o','DisplayName','$Re(\omega_{k_y})$'); hold on;
plot(linear_gr.ky,linear_gr.w_ky(:,end),'--*','DisplayName','$Im(\omega_{k_y})$'); hold on;
plot(linear_gr.ky,linear_gr.ce (:,end),'x','DisplayName','$\epsilon$'); hold on;
+
+ errorbar(linear_gr.ky,linear_gr.avg_g,linear_gr.std_g,...
+ '-o','DisplayName','$Re(\omega_{k_y})$',...
+ 'LineWidth',1.5); hold on;
+ errorbar(linear_gr.ky,linear_gr.avg_w,linear_gr.std_w,...
+ '--*','DisplayName','$Im(\omega_{k_y})$',...
+ 'LineWidth',1.5); hold on;
+
% xlim([min(linear_gr.ky) max(linear_gr.ky)]);
xlabel('$k_y$');
legend('show');
title(DATA.param_title);
- if PLOT > 1
- [YY,XX] = meshgrid(kys,t(its+1:ite));
- figure
- subplot(311)
-% imagesc(t(its:ite),kys,real(w_ky)); set(gca,'YDir','normal');
- pclr = pcolor(XX',YY',real(w_ky)); set(pclr, 'edgecolor','none');
- xlabel('$t$'); ylabel('$k_y$');
- title('Re($\omega$)')
-
- subplot(312)
- pclr = pcolor(XX',YY',imag(w_ky)); set(pclr, 'edgecolor','none');
- xlabel('$t$'); ylabel('$k_y$');
- title('Im($\omega$)')
-
- subplot(313)
- pclr = pcolor(XX',YY',abs(w_ky)); set(pclr, 'edgecolor','none');
- xlabel('$t$'); ylabel('$k_y$');
- title('|$\omega$|')
- end
+
+if PLOT > 1
+ subplot(1,PLOT,2)
+ plot(DATA.Ts3D(its+1:ite),linear_gr.g_ky(Is,:)); hold on;
+ plot(DATA.Ts3D(its+1:ite),linear_gr.w_ky(Is,:));
+ xlabel('t'); ylabel('$\gamma(t),\omega(t)$')
+end
+
+if PLOT > 2
+ subplot(1,PLOT,3)
+ semilogy(DATA.Ts3D,squeeze(abs(DATA.PHI(2,1,DATA.Nz/2,:)))); hold on;
+ xlabel('t'); ylabel('$|\phi_{ky}|(t)$')
end
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index e996e71103944b15f19190d019f12bcb7331d631..b9b99f3ae976c17f8dc5a0735440b5f6c5b9667b 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -8,12 +8,13 @@ for i = 1:numel(dataObjs)
Y_ = [Y_ dataObjs(i).YData];
end
n0 = 1;
+n1 = 26;%numel(X_);
figure;
mvm = @(x) movmean(x,1);
shift = X_(n0);
% shift = 0;
% plot(X_(n0:end),Y_(n0:end));
-plot(mvm(X_(n0:end)-shift),mvm(Y_(n0:end))); hold on;
+plot(mvm(X_(n0:n1)-shift),mvm(Y_(n0:n1))); hold on;
t0 = ceil(numel(X_)*0.5); t1 = numel(X_);
avg= mean(Y_(t0:t1)); dev = std(Y_(t0:t1));
diff --git a/matlab/load/quick_gene_load.m b/matlab/load/quick_gene_load.m
index 0e0e991d9d1e5b8fa44d1c15ac36a67ee6cccc04..3ac095c94f09c8abae00640542a841853413bdd0 100644
--- a/matlab/load/quick_gene_load.m
+++ b/matlab/load/quick_gene_load.m
@@ -41,8 +41,15 @@
% fname ='GENE_LIN_Kn_1.6_KT_0.4_nu_0_32x16.txt';
% fname ='GENE_LIN_Kn_2.5_KT_0.625_nu_0_32x16.txt';
path = '/home/ahoffman/gene/linear_CBC_results/';
-fname = 'CBC_linear.txt';
+% fname = 'CBC_100_20x1x32x32x12_Lv_3_Lw_12.txt';
+% fname = 'CBC_KT_4_20x1x32x32x12_Lv_3_Lw_12.txt';
+% fname = 'CBC_KT_4_20x1x32x64x24_Lv_6_Lw_24.txt';
+% fname = 'CBC_KT_5.3_20x1x32x32x12_Lv_3_Lw_12.txt';
+% fname = 'CBC_KT_5.3_32x1x48x40x16_Lv_3_Lw_12.txt';
+fname = 'CBC_ky_0.3_20x1x32x32x12_Lv_3_Lw_12.txt';
+% fname = 'CBC_ky_0.3_20x1x32x32x12_Lv_3_Lw_12_nuv_1e-3.txt';
data_ = load([path,fname]);
figure
-plot(data_(:,2),data_(:,3),'-dk','DisplayName','GENE');
\ No newline at end of file
+plot(data_(:,2),data_(:,3),'-dk','DisplayName','GENE'); hold on;
+plot(data_(:,2),data_(:,4),'--*k','DisplayName','GENE');
\ No newline at end of file
diff --git a/matlab/plot/spectrum_1D.m b/matlab/plot/spectrum_1D.m
index 14ebd4946581a1faf145706c6eff7bfbf29ce2e5..7288fe1764d4cf00576da25b73dc4f33d587b394 100644
--- a/matlab/plot/spectrum_1D.m
+++ b/matlab/plot/spectrum_1D.m
@@ -159,6 +159,9 @@ else
end
pclr = pcolor(toplot.X,toplot.Y,Gk);
set(pclr, 'edgecolor','none');
+ xlabel('$k_x$');
+ ylabel('$k_y$');
+ title(['GM $|$',fieldname,'$(k_x,k_y)|$ t-averaged']);
end
diff --git a/matlab/setup.m b/matlab/setup.m
index 6228cac6a33725440c569c406db5704100c8259f..db703d0145b228acf16d253ce3ae1f8acfedcc53 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -44,7 +44,7 @@ if MODEL.q_e == 0; SIMID = [SIMID,'_i']; end;
% gradients L_perp/L_x
MODEL.K_n = K_N; % source term kappa for HW
MODEL.K_T = K_T; % Temperature
-MODEL.K_E = K_E; % Electric
+MODEL.K_E = 0; % Electric
MODEL.GradB = GRADB; % Magnetic gradient
MODEL.CurvB = CURVB; % Magnetic curvature
MODEL.lambdaD = LAMBDAD;
@@ -98,7 +98,7 @@ drives_ = [];
if abs(K_N) > 0; drives_ = [drives_,'_kN_',num2str(K_N)]; end;
if abs(K_T) > 0; drives_ = [drives_,'_kT_',num2str(K_T)]; end;
% collision
-coll_ = ['_nu_%0.0e_',CONAME];
+coll_ = ['_nu_%1.1e_',CONAME];
coll_ = sprintf(coll_,NU);
% nonlinear
lin_ = [];
diff --git a/wk/CBC_kT_PJ_scan_quick_run.m b/wk/CBC_kT_PJ_scan_quick_run.m
new file mode 100644
index 0000000000000000000000000000000000000000..4245d7d5502f243fcbd62ea2bf150dfc5d65a54f
--- /dev/null
+++ b/wk/CBC_kT_PJ_scan_quick_run.m
@@ -0,0 +1,47 @@
+KT_a = [3:0.5:7];
+g_max= KT_a*0;
+g_avg= KT_a*0;
+g_std= KT_a*0;
+k_max= KT_a*0;
+
+NU = 0.05;
+DT = 2e-3;
+TMAX = 50;
+ky_ = 0.3;
+SIMID = 'linear_CBC_KT_scan_ky=0.3_CLOS_0_DGGK_0.05'; % Name of the simulation
+% SIMID = 'linear_CBC_PJ_scan_KT_4_ky=0.3_CLOS_0_TMAX_100'; % Name of the simulation
+RUN = 1;
+figure
+
+for P = [20]
+
+i=1;
+for K_T = KT_a
+
+ quick_run
+
+ g_max(i) = gmax;
+ k_max(i) = kmax;
+
+ g_avg(i) = lg.avg_g;
+ g_std(i) = lg.std_g;
+
+ i = i + 1;
+end
+%%
+
+% plot(KT_a,max(g_max,0));
+y_ = g_avg;
+e_ = g_std;
+
+y_ = y_.*(y_-e_>0);
+e_ = e_ .* (y_>0);
+errorbar(KT_a,y_,e_,...
+ 'LineWidth',1.2,...
+ 'DisplayName',['(',num2str(P),',',num2str(P/2),')']);
+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/CBC_nu_CO_scan_quick_run.m b/wk/CBC_nu_CO_scan_quick_run.m
new file mode 100644
index 0000000000000000000000000000000000000000..6b9bc4fa9da0ccc5147004ef0518f97f5996685b
--- /dev/null
+++ b/wk/CBC_nu_CO_scan_quick_run.m
@@ -0,0 +1,48 @@
+% NU_a = [0.05 0.15 0.25 0.35 0.45];
+NU_a = [0:0.05:0.5];
+g_max= NU_a*0;
+g_avg= NU_a*0;
+g_std= NU_a*0;
+k_max= NU_a*0;
+CO = 'LR';
+
+K_T = 5.3;
+DT = 2e-3;
+TMAX = 30;
+ky_ = 0.3;
+SIMID = 'linear_CBC_nu_scan_ky=0.3_CLOS_0_LRGK'; % Name of the simulation
+RUN = 1;
+figure
+
+for P = [2 4 6 10 12 20]
+
+i=1;
+for NU = NU_a
+
+ quick_run
+
+ g_max(i) = gmax;
+ k_max(i) = kmax;
+
+ g_avg(i) = lg.avg_g;
+ g_std(i) = lg.std_g;
+
+ i = i + 1;
+end
+%%
+
+% plot(KT_a,max(g_max,0));
+y_ = g_avg;
+e_ = g_std;
+
+y_ = y_.*(y_-e_>0);
+e_ = e_ .* (y_>0);
+errorbar(NU_a,y_,e_,...
+ 'LineWidth',1.2,...
+ 'DisplayName',['(',num2str(P),',',num2str(P/2),')']);
+hold on;
+title(['$\kappa_T=$',num2str(K_T),' $k_y=$',num2str(ky_),' (CLOS = 0)']);
+legend('show'); xlabel('$\nu_{DGGK}$'); ylabel('$\gamma$');
+drawnow
+end
+
diff --git a/wk/Dimits_2000_fig3.m b/wk/Dimits_2000_fig3.m
new file mode 100644
index 0000000000000000000000000000000000000000..96c8b4b25f0bab29082531886a86f23029bc0612
--- /dev/null
+++ b/wk/Dimits_2000_fig3.m
@@ -0,0 +1,28 @@
+%192x96x16x3x2 simulation
+kT_Xi_GM_32 = ...
+ [7.0 3.6;...
+ 6.0 2.6;...
+ 5.0 1.1;...
+ ];
+%128x64x16x5x3 simulation
+kT_Xi_GM_53 = ...
+ [7.0 4.6;...
+ 5.3 1.9;...
+ 5.0 1.5;...
+ ];
+%128x64x16x24x12 simulation
+kT_Xi_GN = ...
+ [7.0 5.0;...
+ 5.3 2.4;...
+ 4.5 nan;...
+ ];
+
+figure;
+plot(kT_Xi_GM_32(:,1), kT_Xi_GM_32(:,2),'o-','DisplayName','192x96x16x3x2'); hold on
+plot(kT_Xi_GM_53(:,1), kT_Xi_GM_53(:,2),'o-','DisplayName','128x64x16x5x3'); hold on
+plot(kT_Xi_GN(:,1), kT_Xi_GN(:,2),'o--k','DisplayName','GENE 128x64x16x24x12');
+xlabel('$\kappa_T$'); ylabel('$\chi_i$');
+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
new file mode 100644
index 0000000000000000000000000000000000000000..9671fc0fe8d54e9f14220bca4fb29cebd40709a4
--- /dev/null
+++ b/wk/ITG_peak_KT_53_k_035_scan.m
@@ -0,0 +1,6 @@
+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/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index 7721a364a50f299fed8bed4af9a168d08747bdac..1ba1f5aafbb5318c25dc6bba07ef9feb9d10322b 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -11,7 +11,7 @@ system(['mkdir -p ',LOCALDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
system(CMD);
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 20;
+JOBNUMMIN = 00; JOBNUMMAX = 10;
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
data.localdir = LOCALDIR;
data.FIGDIR = LOCALDIR;
@@ -23,7 +23,7 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-options.TAVG_0 = 0.7*data.Ts3D(end); data.scale = 1;%(1/data.Nx/data.Ny)^2;
+options.TAVG_0 = 0.7*data.Ts3D(end); data.scale = 1;%/(data.Nx*data.Ny)^2;
options.TAVG_1 = data.Ts3D(end); % Averaging times duration
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)
@@ -55,8 +55,8 @@ options.PLAN = 'xz';
% options.PLAN = 'sx';
options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
-options.TIME = data.Ts3D(1:2:end);
-% options.TIME = [750:1:1000];
+% options.TIME = data.Ts3D(1:10:end);
+options.TIME = [200:2:400];
data.EPS = 0.1;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
@@ -67,18 +67,18 @@ if 0
% Options
options.INTERP = 0;
options.POLARPLOT = 0;
-options.AXISEQUAL = 1;
+options.AXISEQUAL = 0;
% options.NAME = '\phi';
% options.NAME = 'n_e';
options.NAME = 'N_i^{00}';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'kxz';
+options.PLAN = 'kxky';
% options.NAME 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [120 150 155];
+options.TIME = [20 60 200 400 500];
data.a = data.EPS * 2e3;
fig = photomaton(data,options);
% save_figure(data,fig)
@@ -98,12 +98,12 @@ end
if 0
%% Kinetic distribution function sqrt(<f_a^2>xy) (GENE vsp)
-options.SPAR = linspace(-3,3,32)+(6/127/2);
-options.XPERP = linspace( 0,6,32);
-% options.SPAR = gene_data.vp';
-% options.XPERP = gene_data.mu';
-options.iz = 'avg';
-options.T = [600];
+% options.SPAR = linspace(-3,3,32)+(6/127/2);
+% options.XPERP = linspace( 0,6,32);
+options.SPAR = gene_data.vp';
+options.XPERP = gene_data.mu';
+options.iz = 1;
+options.T = [500 1000];
options.PLT_FCT = 'contour';
options.ONED = 0;
options.non_adiab = 0;
@@ -119,7 +119,7 @@ if 0
options.P2J = 0;
options.ST = 1;
options.PLOT_TYPE = 'space-time';
-options.NORMALIZED = 1;
+options.NORMALIZED = 0;
options.JOBNUM = 0;
options.TIME = [1000];
options.specie = 'i';
@@ -131,11 +131,11 @@ end
if 0
%% Time averaged spectrum
-options.TIME = [300 600];
+options.TIME = [100 500];
options.NORM =1;
% options.NAME = '\phi';
-% options.NAME = 'N_i^{00}';
-options.NAME ='\Gamma_x';
+options.NAME = 'N_i^{00}';
+% options.NAME ='\Gamma_x';
options.PLAN = 'kxky';
options.COMPZ = 'avg';
options.OK = 0;
@@ -169,7 +169,7 @@ options.NORMALIZED = 0;
options.K2PLOT = 1;
options.TIME = [0:160];
options.NMA = 1;
-options.NMODES = 15;
+options.NMODES = 1;
options.iz = 'avg';
fig = mode_growth_meter(data,options);
save_figure(data,fig,'.png')
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 4b20e0b514f73158bf4046250ca94836175bf6de..57ba009efb1863e599a7171de4b6abdd7680cd01 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -13,6 +13,9 @@
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
% folder = '/misc/gene_results/ZP_kn_2.5_large_box/';
folder = '/misc/gene_results/CBC/128x64x16x24x12/';
+% folder = '/misc/gene_results/CBC/196x96x20x32x16_00/';
+% folder = '/misc/gene_results/CBC/128x64x16x6x4/';
+% folder = '/misc/gene_results/CBC/KT_5.3_128x64x16x24x12_00/';
gene_data = load_gene_data(folder);
gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1
@@ -101,7 +104,7 @@ end
if 0
%% Show f_i(vpar,mu)
-options.times = 200:300;
+options.times = 20:80;
options.specie = 'i';
options.PLT_FCT = 'contour';
options.folder = folder;
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 9b6139f7c5a8e86ee113fe1187204fb31c77ee5d..45ec1e032f649ed9860c25e02d54ef67e9f1001e 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -31,7 +31,7 @@ helazdir = '/home/ahoffman/HeLaZ/';
% outfile ='Zpinch_rerun/Kn_2.5_256x128x5x3';
% outfile ='Zpinch_rerun/Kn_2.5_312x196x5x3_Lx_400_Ly_200';
% outfile ='Zpinch_rerun/Kn_2.5_256x64x5x3';
-outfile ='Zpinch_rerun/Kn_2.0_200x48x9x5_large_box';
+% outfile ='Zpinch_rerun/Kn_2.0_200x48x9x5_large_box';
% outfile ='Zpinch_rerun/Kn_2.0_256x64x9x5_Lx_240_Ly_120';
% outfile ='Zpinch_rerun/Kn_1.6_256x128x7x4';
% outfile ='Zpinch_rerun/Kn_1.6_200x48x11x6';
@@ -41,5 +41,12 @@ outfile ='Zpinch_rerun/Kn_2.0_200x48x9x5_large_box';
% outfile = 'CBC/64x32x16x5x3';
% outfile = 'CBC/64x128x16x5x3';
% outfile = 'CBC/128x64x16x5x3';
-outfile = 'CBC/64x64x16x3x2';
+% outfile = 'CBC/192x96x16x3x2';
+% outfile = 'CBC/128x64x16x5x3';
+
+outfile = 'CBC/kT_scan_128x64x16x5x3';
+% outfile = 'CBC/kT_scan_192x96x16x3x2';
+
+%% 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';
run analysis_HeLaZ
diff --git a/wk/quick_run.m b/wk/quick_run.m
index e0bcee2e6129c604857f2977ce59687d6a8c52a2..8a7931f43108c194a9f2d02f1057d287d1b686cc 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -3,35 +3,39 @@
% from matlab framework. It is meant to run only small problems in linear
% for benchmark and debugging purpose since it makes matlab "busy"
%
-RUN = 1; % To run or just to load
+% RUN = 1; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
HELAZDIR = '/home/ahoffman/HeLaZ/';
EXECNAME = 'helaz3';
-% EXECNAME = 'helaz3_shear'; %verified version
+% EXECNAME = 'helaz3_dbg';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.01; % Collision frequency
+% NU = 0.0; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_N = 2.22;%2.0; % Density gradient drive
-K_T = 6.92;%0.25*K_N; % Temperature '''
-K_E = 0.0; % Electrostat '''
+K_N = 2.22;%2.0; % ion Density gradient drive
+K_Ne = K_N; % ele Density gradient drive
+% K_T = 4.0;%0.25*K_N; % Temperature '''
+K_Te = 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 = 0e-1; % electron plasma beta
%% GRID PARAMETERS
-PMAXE = 4; % Hermite basis size of electrons
-JMAXE = 2; % Laguerre "
-PMAXI = 4; % " ions
-JMAXI = 2; % "
+% P = 12;
+J = P/2;
+PMAXE = P; % Hermite basis size of electrons
+JMAXE = J; % Laguerre "
+PMAXI = P; % " ions
+JMAXI = J; % "
NX = 12; % real space x-gridpoints
-NY = 8; % '' y-gridpoints
+NY = 2; % '' y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain
-LY = 2*pi/0.1; % Size of the squared frequency domain
+% LY = 2*pi/0.3; % Size of the squared frequency domain
+LY = 2*pi/ky_;
NZ = 16; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
@@ -42,8 +46,8 @@ Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear (Not implemented yet)
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 20; % Maximal time unit
-DT = 5e-3; % Time step
+% TMAX = 100; % Maximal time unit
+% DT = 2e-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
@@ -51,15 +55,16 @@ SPS5D = 1/5; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
-SIMID = 'linear_CBC'; % Name of the simulation
+% SIMID = 'linear_CBC'; % Name of the simulation
+% SIMID = 'scan'; % Name of the simulation
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
+% CO = 'SG';
+GKCO = 1; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
-CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))s
+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
@@ -96,7 +101,7 @@ 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,' 2 1 3 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
end
@@ -112,7 +117,7 @@ data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from
if 1
%% linear growth rate (adapted for 2D zpinch and fluxtube)
trange = [0.5 1]*data.Ts3D(end);
-nplots = 1;
+nplots = 0;
lg = compute_fluxtube_growth_rate(data,trange,nplots);
[gmax, kmax] = max(lg.g_ky(:,end));
[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
@@ -157,17 +162,16 @@ options.kzky = 0;
[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
save_figure(data,fig)
end
-
if 0
%% Mode evolution
-options.NORMALIZED = 1;
+options.NORMALIZED = 0;
options.K2PLOT = 1;
-options.TIME = [0 1]*data.Ts3D(end);
+options.TIME = [0:1000];
options.NMA = 1;
options.NMODES = 1;
-options.iz = 9;
+options.iz = 'avg';
fig = mode_growth_meter(data,options);
-save_figure(gbms_dat,fig)
+save_figure(data,fig,'.png')
end
diff --git a/wk/scan_quick_run.m b/wk/scan_quick_run.m
new file mode 100644
index 0000000000000000000000000000000000000000..4245d7d5502f243fcbd62ea2bf150dfc5d65a54f
--- /dev/null
+++ b/wk/scan_quick_run.m
@@ -0,0 +1,47 @@
+KT_a = [3:0.5:7];
+g_max= KT_a*0;
+g_avg= KT_a*0;
+g_std= KT_a*0;
+k_max= KT_a*0;
+
+NU = 0.05;
+DT = 2e-3;
+TMAX = 50;
+ky_ = 0.3;
+SIMID = 'linear_CBC_KT_scan_ky=0.3_CLOS_0_DGGK_0.05'; % Name of the simulation
+% SIMID = 'linear_CBC_PJ_scan_KT_4_ky=0.3_CLOS_0_TMAX_100'; % Name of the simulation
+RUN = 1;
+figure
+
+for P = [20]
+
+i=1;
+for K_T = KT_a
+
+ quick_run
+
+ g_max(i) = gmax;
+ k_max(i) = kmax;
+
+ g_avg(i) = lg.avg_g;
+ g_std(i) = lg.std_g;
+
+ i = i + 1;
+end
+%%
+
+% plot(KT_a,max(g_max,0));
+y_ = g_avg;
+e_ = g_std;
+
+y_ = y_.*(y_-e_>0);
+e_ = e_ .* (y_>0);
+errorbar(KT_a,y_,e_,...
+ 'LineWidth',1.2,...
+ 'DisplayName',['(',num2str(P),',',num2str(P/2),')']);
+hold on;
+title(['Linear CBC $K_T$ threshold $k_y=$',num2str(ky_),' (CLOS = 1)']);
+legend('show'); xlabel('$K_T$'); ylabel('$\gamma$');
+drawnow
+end
+