diff --git a/wk/Ch7_HF_analysis.m b/wk/Ch7_HF_analysis.m
new file mode 100644
index 0000000000000000000000000000000000000000..750ec4b91f6f3a950b95b72ff7f7e2d61ba75ca1
--- /dev/null
+++ b/wk/Ch7_HF_analysis.m
@@ -0,0 +1,45 @@
+PARTITION = '/misc/gyacomo23_outputs/paper_3/';
+switch 3
+case 1
+ SIM_SET_NAME = 'Multi-scale';
+ E_FLUX = 1;
+ FOLDER = 'DIIID_fullphys_rho95_geom_scan/multi_scale/multi_scale_3x2x768x192x24';
+ % FOLDER = 'DIIID_fullphys_rho95_geom_scan/multi_scale/multi_scale_5x2x768x192x24';
+case 2
+ SIM_SET_NAME = 'Ion-scale';
+ E_FLUX = 1;
+ FOLDER = 'DIIID_fullphys_rho95_geom_scan/ion_scale/ion_scale_5x2x256x64x32_tau_1_RN_0';
+case 3
+ SIM_SET_NAME = 'Adiab. e.';
+ E_FLUX = 0;
+ FOLDER = 'DIIID_adiab_e_rho95_geom_scan/5x2x256x64x32_tau_1_RN_0/';
+case 4
+ SIM_SET_NAME = 'Adiab. e.';
+ E_FLUX = 0;
+ FOLDER = 'DIIID_adiab_e_rho95_geom_scan/5x2x256x64x32_tau_1_RN_0/';
+end
+
+GEOM = {'NT','0T','PT'};
+COLO = {'b','k','r'};
+figure
+for i = 1:3
+ OPTION = GEOM{i};
+ datadir = [PARTITION,FOLDER,'/',OPTION];
+ out = read_flux_out_XX(datadir,0);
+
+ [ts, is] = sort(out.t);
+ Pxis = out.Pxi(is);
+ Qxis = out.Qxi(is);
+
+ if E_FLUX
+ Pxes = out.Pxe(is);
+ Qxes = out.Qxe(is);
+ subplot(211)
+ plot(ts,Qxes,COLO{i}); hold on;
+ subplot(212)
+ title(SIM_SET_NAME)
+ end
+ plot(ts,Qxis,COLO{i}); hold on;
+end
+title(SIM_SET_NAME)
+
diff --git a/wk/HEL_lin_solver.m b/wk/HEL_GM_eig_solver.m
similarity index 82%
rename from wk/HEL_lin_solver.m
rename to wk/HEL_GM_eig_solver.m
index 752fb87f5e8c45cd9aff50f6f19cc2df5f8a8f92..457f4417847a408aead1211f6273028cd603cc04 100644
--- a/wk/HEL_lin_solver.m
+++ b/wk/HEL_GM_eig_solver.m
@@ -1,24 +1,33 @@
% We formulate the linear system as dN/dt + A*N = 0
+% parameters
+q = 1; % Safety factor
+Jxyz = 1; % Jacobian
+hatB = 1; % normalized B field
+dzlnB = 0; % B parallel derivative
+ddz = 0; % parallel derivative operator
+tau = 0.001; % ion-electron temperature ratio
+IVAN = 0; % flag for Ivanov scaling
+kT = 7; % scaled temperature gradient
+chi = 0.16; % scaled collisionality
+kx_a = 0;linspace(-2.0,2.0,256); % radial fourier grid
+ky_a = linspace(0.01,5,256); % binormal fourier grid
+
+% translate into parameters in GM formulation
+RN = 0;
+RT = kT*2*q/tau;
+NU = chi*3/2/tau/(4-IVAN*2.25);
+MU = 0.0;
-q = 1;
-Jxyz = 1; hatB = 1; dzlnB = 0;
-ddz = 0;
-tau = 0.001;
-IVAN = 1;
%ordering
O1_n = 1;
O1_u = 1-IVAN;
O1_Tpar = 1-IVAN;
O1_Tper = 1-IVAN;
-RN = 0;
-RT = 1*2*q/tau;
-NU = 0.1*3/2/tau/(4-IVAN*2.25);
-MU = 0.0;
-
-kx_a = 0;linspace(-2.0,2.0,256);
-ky_a = linspace(0.01,3.5,256);
+% array of most unstable growth rates and corresponding frequencies
g_a = zeros(numel(kx_a),numel(ky_a)); w_a = g_a;
+
+% Evaluation of the matrices and solver
for i = 1:numel(kx_a)
for j = 1:numel(ky_a)
kx = kx_a(i);
@@ -116,26 +125,8 @@ for j = 1:numel(ky_a)
CDG(4,3) =-2/3*sqrt(2)*K0*(K0-2*K1*O1_Tper);
CDG(4,4) = 4/3*(K0-2*K1*O1_Tper)^2 + 2*tau*(K0-K1)^2*lperp*O1_Tper;
CDG(4,5) =-2/3*q/tau*(K0-K1)*(3*tau*K0^2*lperp-K0*K1*(4+3*tau*lperp)+K1^2*(8+3*tau*lperp));
- if 1
- C = NU*(CLB + CDG);
- else %to check
- C = zeros(4,5);
- C(1,1) = -2/3*tau*lperp*(4*tau*lperp);
- C(1,3) = -2/3*tau*lperp*(sqrt(2)-2*sqrt(2)*tau*lperp);
- C(1,4) = -2/3*tau*lperp*(1-tau*lperp);
- C(1,5) = -2/3*tau*lperp*(5*q*lperp - 11*q*tau*lperp^2);
- C(2,2) = -(1+2*tau*lperp);
- C(3,1) = -2/3*(sqrt(2)*tau*lperp);
- C(3,3) = -2/3*(2+5*tau*lperp);
- C(3,4) = -2/3*(sqrt(2)-4*sqrt(2)*tau*lperp);
- C(3,5) = -2/3*(2*sqrt(2)*q*lperp+8*sqrt(2)*q*tau*lperp^2);
- C(4,1) = -2/3*(tau*lperp - tau^2*lperp^2);
- C(4,3) = -2/3*(sqrt(2)*4*sqrt(2)*tau*lperp);
- C(4,4) = -2/3*(1+12*tau*lperp);
- C(4,5) = -2/3*(2*q*lperp+9*q*tau*lperp^2);
- C = NU*C;
+ C = NU*(CLB + CDG);
- end
% Numerical diffusion
Diff = MU*lperp^2*eye(4);
@@ -166,7 +157,11 @@ for j = 1:numel(ky_a)
gamma =-real(lambda);
omega = imag(lambda);
[g_a(i,j),idx] = max(gamma);
- % w_a(i,j) = omega(idx);
+ % gsorted = sort(gamma);
+ % g_a(i,j) = gsorted(1);
+ w_a(i,j) = omega(idx);
+ % wsorted = sort(omega);
+ % g_a(i,j) = wsorted(1);
w_a(i,j) = max(omega);
end
end
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index 215a920d5753933fa524a6cdb26377480ee81713..4e0e5814654e347f9d832e3af220307ef7567f15 100644
--- a/wk/fast_analysis.m
+++ b/wk/fast_analysis.m
@@ -7,7 +7,7 @@ default_plots_options
% Partition of the computer where the data have to be searched
% PARTITION='/Users/ahoffmann/gyacomo/results/paper_3/';
PARTITION='/misc/gyacomo23_outputs/paper_3/';
-% PARTITION = '';
+% PARTITION = '../results/paper_3/';
%% Paper 3
% resdir = 'DTT_rho85/3x2x192x48x32';
% resdir = 'DTT_rho85/3x2x192x48x32_NT';
@@ -17,88 +17,168 @@ PARTITION='/misc/gyacomo23_outputs/paper_3/';
% resdir = 'LM_DIIID_rho95/3x2x512x92x32';
% resdir = 'DIIID_LM_rho90/3x2x256x128x32';
% resdir = 'DTT_rho85_geom_scan/P8_J4_delta_nuDGGK_conv_test/delta_-0.3_nu_0.9';
-resdir = 'NT_DIIID_Austin2019_rho95/3x2x256x64x32';
-% resdir = '../testcases/cyclone_example';
+% resdir = 'NT_DIIID_Austin2019_rho95/3x2x256x64x32';
+
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/largerbox_moremodes';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/PT/3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/NT/3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/VNT/3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/VPT/3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/CIRC/3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/ELONG/3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/PT/lin_3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/NT/lin_3x2x128x32x32';
+
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/huge';
+
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_1000/PT/3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_1000/NT/3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_1000/PT/lin_3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_1000/NT/lin_3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_750/PT/lin_3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_750/NT/lin_3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_250/PT/lin_3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_250/NT/lin_3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-3/huge';
+
+% resdir = 'DIIID_rho95_cold_ions_tau1e-2/PT/3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e-2/NT/3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e0/PT/3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e0/PT/5x3x192x48x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e0/NT/3x2x128x32x32';
+% resdir = 'DIIID_rho95_cold_ions_tau1e0/NT/5x3x192x48x32';
+% resdir = 'DTT_rho85_cold_ions_tau1e-3/NT/3x2x128x32x32';
+% % resdir = '../testcases/cyclone_example';
% resdir = '../testcases/CBC_ExBshear/std';
% resdir = '../results/paper_3/HM_DTT_rho98/3x2x128x64x64';
- %%
-J0 = 00; J1 = 10;
-% Load basic info (grids and time traces)
+% resdir = 'DIIID_cold_ions_rho95_geom_scan/3x2x192x48x32_RT_1000_eps_q0_scan/NT/eps_0.35_q0_4.0';
+% resdir = 'DTT_rho85_geom_scan/P2_J1_PT_sfact_shear_scan/shear_2.7_q0_-2.9';
+% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_kT_3500/128x32x24';
+% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_kT_3500/256x64x24';
+% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_colless';
+% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_kT_2000/lin_128x32x24';
+% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_kT_2000/128x32x24';
+% PARTITION = ''; resdir ='../results/HEL_CBC/CBC_21/128x32x24';
+% PARTITION = ''; resdir ='../results/HEL_CBC/192x48x24';
+% PARTITION = ''; resdir ='../testcases/Ivanov_2020';
+% PARTITION = ''; resdir ='../testcases/Hasegawa_Wakatani';
+% resdir = 'HEL_CBC/256x92x24_max_trunc';
+
+% resdir ='HEL_CBC/192x48x24';
+% resdir ='HEL_CBC/256x92x24'3;
+% resdir ='HEL_CBC/256x256x32/k_N__0.0_k_T__1750';
+
+PARTITION = '/misc/gyacomo23_outputs/paper_3/DIIID_rho_95_Tstudy/';
+% resdir = 'multi_scale_3x2x512x128x24';
+% resdir = 'multi_scale_3x2x512x128x24_larger_box';
+% resdir = 'multi_scale_3x2x768x192x24/continue_with_gradN_and_tau';
+% resdir = 'multi_scale/multi_scale_5x2x768x192x24/PT';
+% resdir = 'NT';
+% resdir = 'electron_scale_3x2x256x64x24';
+% resdir = 'electron_scale_3x2x128x64x24';
+% resdir = 'ion_scale_3x2x192x48x32_larger_box';
+% resdir = 'ion_scale_3x2x256x64x32'; % PT and NT also
+% resdir = 'ion_scale/ion_scale_5x2x256x64x32_tau_1_RN_0/NT';
+% resdir = 'adiab_e/5x2x256x64x32_tau_1_RN_0/NT';
+
+resdir = 'multi_scale/multi_scale_3x2x768x192x24/NT';
+% resdir = 'multi_scale/multi_scale_3x2x512x128x24_larger_box/NT';
+% resdir = 'multi_scale/multi_scale_3x2x768x192x24/continue_with_gradN_and_tau';
+% resdir = 'ion_scale/ion_scale_5x2x256x64x32_tau_1_RN_0/NT';
+% resdir = 'adiab_e/5x2x256x64x32_tau_1_RN_0/0T';
+% resdir = 'hot_electrons/hot_electrons_256x64x24/0T';
+
+
+
+% PARTITION = '/misc/gyacomo23_outputs/paper_3/';
+% resdir = 'DIIID_HEL_rho95/PT';
+
DATADIR = [PARTITION,resdir,'/'];
+%%
+J0 = 00; J1 = 20;
+
+% Load basic info (grids and time traces)
data = {};
data = compile_results_low_mem(data,DATADIR,J0,J1);
-
+[data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00');
+data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+try
+data.Ne00 = reshape(data.Na00(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+catch
+end
+[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
+if 1
+ %%
+[data.TEMP, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'temp');
+% [data.UPAR, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'upar');
+[data.DENS, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'dens');
+data.TEMP_I = reshape(data.TEMP(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+% data.UPAR_I = reshape(data.UPAR(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+data.DENS_I = reshape(data.DENS(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+if data.inputs.Na > 1
+ data.TEMP_E = reshape(data.TEMP(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+ data.DENS_E = reshape(data.DENS(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+ data.Ne00 = reshape(data.Na00(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+end
+end
if 1
%% Plot transport and phi radial profile
-[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
+% [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
% [data.PSI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'psi');
-options.TAVG_0 = 100;
-options.TAVG_1 = 1000;
+options.TAVG_0 = data.Ts3D(end)/2;
+options.TAVG_1 = data.Ts3D(end);
options.NCUT = 5; % 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)
+% options.ST_FIELD = 'n_e'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+% options.ST_FIELD = 'u_i'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+% options.ST_FIELD = 'T_i'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+% options.ST_FIELD = 'n_i T_i'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+% options.ST_FIELD = 'Q_{xi}'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+% options.ST_FIELD = 'G_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+% options.ST_FIELD = 'w_{Ez}'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+% options.ST_FIELD = 'v_{Ey}'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+% options.ST_FIELD = 'N_i^{00}'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
options.INTERP = 0;
options.RESOLUTION = 256;
plot_radial_transport_and_spacetime(data,options);
end
-if 0
-%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% Options
-[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
-[data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00');
-data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
-options.INTERP = 1;
-options.POLARPLOT = 0;
-options.BWR = 0; % bluewhitered plot or gray
-options.CLIMAUTO = 1; % adjust the colormap auto
-% options.NAME = '\phi';
-% options.NAME = '\phi^{NZ}';
-% options.NAME = '\omega_z';
-options.NAME = 'N_i^{00}';
-% options.NAME = 's_{Ey}';
-% options.NAME = 'n_i^{NZ}';
-% options.NAME = 'Q_x';
-% options.NAME = 'n_i';
-% options.NAME = 'n_i-n_e';
-options.PLAN = 'kxky';
-% options.NAME = 'f_i';
-% options.PLAN = 'sx';
-options.COMP = 'avg';
-% options.TIME = data.Ts5D(end-30:end);
-options.TIME = data.Ts3D(1:1:end);
-% options.TIME = [0:1500];
-data.EPS = 0.1;
-data.a = data.EPS * 2000;
-options.RESOLUTION = 256;
-create_film(data,options,'.gif')
-end
-
-if 0
-%% field snapshots
+if 1
+%% 2D field snapshots
% Options
-[data.Na00, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'Na00');
-[data.PHI, data.Ts3D] = compile_results_3D(data.folder,J0,J1,'phi');
-data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
-
options.INTERP = 0;
options.POLARPLOT = 0;
options.AXISEQUAL = 0;
options.NORMALIZE = 0;
options.LOGSCALE = 0;
options.CLIMAUTO = 1;
-options.NAME = 'N_i^{00}';
-% options.NAME = 's_{Ey}';
+options.TAVG = 1;
+% options.NAME = ['N_e^{00}'];
+% options.NAME = 'n_e';
+% % options.NAME = 'u_i';
+% options.NAME = 'T_i';
+% options.NAME = 'Q_{xe}';
+options.NAME = 'v_{Ey}';
+% options.NAME = 'w_{Ez}';
+% options.NAME = '\omega_z';
% options.NAME = '\phi';
-options.PLAN = 'yz';
-options.COMP = 'avg';
-options.TIME = [50 100];
-% options.TIME = data.Ts3D(1:2:end);
+% options.NAME = 'n_i-n_e';
+loc =11;
+[~,i_] = min(abs(loc - data.grids.y));
+options.COMP =i_;
+% options.PLAN = '3D';
+options.PLAN = 'xy'; options.COMP =floor(data.grids.Nz/2)+1;
+% options.PLAN = 'yz'; options.COMP ='avg';
+% options.COMP ='avg';
+options.XYZ =[-11 20 0];
+options.TIME = [90:150];
options.RESOLUTION = 256;
fig = photomaton(data,options);
-colormap(gray)
+% colormap(gray)
clim('auto')
% save_figure(data,fig)
end
@@ -107,42 +187,147 @@ if 0
profiler(data)
end
-if 0
+if 1
+%% Mode evolution
+% [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
+% [data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00');
+% data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+% data.Ne00 = reshape(data.Na00(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+
+options.NORMALIZED = 0;
+options.TIME = data.Ts3D;
+options.KX_TW = [0.1 2.5]; %kx Growth rate time window
+options.KY_TW = [0.1 2.5]; %ky Growth rate time window
+options.NMA = 1;
+options.NMODES = 64;
+options.iz = 'avg'; % avg or index
+options.ik = 1; % sum, max or index
+options.fftz.flag = 0;
+options.FIELD = 'Ni00';
+% options.FIELD = 'phi';
+% options.FIELD = 'T_i';
+options.GOK2 = 0;
+options.SHOWFIG = 1;
+[fig, wkykx, ekykx] = mode_growth_meter(data,options);
+%%
+kx = (1:data.grids.Nx/2)'*2*pi/data.fort_00.GRID.Lx;
+ky = (1:data.grids.Ny/2)'*2*pi/data.fort_00.GRID.Ly;
+gkxky = real(wkykx(2:end,1:data.grids.Nx/2))';
+gkxky(isnan(gkxky)) =0;
+gkxky(isinf(gkxky)) =0;
+% gkxky(gkxky<0) =0;
+% gkxky = imgaussfilt(gkxky,1);
+%
+wkxky = imag(wkykx(2:end,1:data.grids.Nx/2))';
+wkxky(isnan(wkxky)) =0;
+wkxky(isinf(wkxky)) =0;
+% wkxky(wkxky<0) =0;
+% wkxky = imgaussfilt(wkxky,1.5);
+%
+figure;
+subplot(121)
+ contourf(kx,ky,gkxky',10)
+ % clim(0.5*[0 1]);
+ % colormap(bluewhitered); colorbar;
+ xlim([0.025 1]);
+ xlabel('$k_x\rho_s$'); ylabel('$k_y\rho_s$')
+subplot(122)
+ contourf(kx,ky,wkxky',10)
+ % clim(1*[0 1]);
+ % colormap(bluewhitered); colorbar
+ xlim([0.025 1]);
+ xlabel('$k_x\rho_s$'); ylabel('$k_y\rho_s$')
+% save_figure(data,fig,'.png')
+end
+
+if 1
%% Hermite-Laguerre spectrum
-[data.Napjz, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Napjz');
+[data.Napjz, data.Ts3D] = compile_results_3Da(DATADIR,0,10,'Napjz');
% [data.Napjz, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'Nipjz');
-options.ST = 0;
+options.ST = 1;
options.NORMALIZED = 0;
-options.LOGSCALE = 1;
+options.LOGSCALE = 0;
options.FILTER = 0; %filter the 50% time-average of the spectrum from
options.TAVG_2D = 0; %Show a 2D plot of the modes, 50% time averaged
options.TAVG_2D_CTR= 0; %make it contour plot
+options.TWINDOW = [6 20];
fig = show_moments_spectrum(data,options);
end
-if 0
-%% Mode evolution
-[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
-[data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00');
-data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+if (0 && NZ>4)
+%% Ballooning plot
+% [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
+if data.inputs.BETA > 0
+[data.PSI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'psi');
+end
+options.time_2_plot = [25 100];
+options.kymodes = 1.5;
+options.normalized = 1;
+options.PLOT_KP = 0;
+% options.field = 'phi';
+options.SHOWFIG = 1;
+[fig, chi, phib, psib, ~] = plot_ballooning(data,options);
+end
-options.NORMALIZED = 0;
-options.TIME = data.Ts3D;
-options.KX_TW = [ 0 20]; %kx Growth rate time window
-options.KY_TW = [ 0 50]; %ky Growth rate time window
-options.NMA = 1;
-options.NMODES = 50;
-options.iz = 'avg'; % avg or index
-options.ik = 9; % sum, max or index
-options.fftz.flag = 0;
-% options.FIELD = 'Ni00';
-options.FIELD = 'phi';
-options.GOK2 = 0;
-fig = mode_growth_meter(data,options);
-% save_figure(data,fig,'.png')
+if 1
+%% 1D spectral plot
+options.TIME = [30 80]; % averaging time window
+% options.NAME = ['N_i^{00}'];
+% options.NAME = 'n_i';
+% options.NAME = 'T_i';
+% options.NAME = 'Q_{xi}';
+% options.NAME = 's_{Ey}';
+options.NAME = '\psi';
+options.NORMALIZE = 0;
+[fig] = plot_spectrum(data,options);
end
+if 0
+%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Options
+% [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
+% [data.PSI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'psi');
+% [data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00');
+% data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+% data.Ne00 = reshape(data.Na00(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+% [data.DENS, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'dens');
+% data.DENS_I = reshape(data.DENS(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+% [data.TEMP, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'temp');
+% data.TEMP_I = reshape(data.TEMP(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+options.INTERP = 0;
+options.POLARPLOT = 0;
+options.BWR = 1; % bluewhitered plot or gray
+options.CLIMAUTO = 0; % adjust the colormap auto
+% options.NAME = '\phi';
+% options.NAME = 'w_{Ez}';
+% options.NAME = '\psi';
+options.NAME = 'T_i';
+% options.NAME = '\phi^{NZ}';
+% options.NAME = ['N_e^{00}'];
+% options.NAME = ['N_i^{00}'];
+options.PLAN = 'xy';
+% options.PLAN = '3D';
+% options.XYZ =[-11 20 0];
+% options.COMP = 'avg';
+options.COMP = floor(data.grids.Nz/2+1);
+options.TIME = data.Ts3D(1:1:end);
+% options.TIME = [0:1500];
+data.EPS = 0.1;
+data.a = data.EPS * 2000;
+options.RESOLUTION = 256;
+options.FPS = 12;
+options.RMAXIS = 0;
+create_film(data,options,'.gif');
+end
+
+if 0
+%% Metric infos
+options.SHOW_FLUXSURF = 1;
+options.SHOW_METRICS = 1;
+[fig, geo_arrays] = plot_metric(data,options);
+end
+
if 0
%% Study singular values
[data.SV_ky_pj, data.Ts2D] = compile_results_2D(DATADIR,J0,J1,'sv_ky_pj');
diff --git a/wk/lin_run_script.m b/wk/lin_run_script.m
index f872fc8767ddeeb77a2475406efafb4050956861..7e0eb7271fab42b84a51e3fea2c12c986b5eee4b 100644
--- a/wk/lin_run_script.m
+++ b/wk/lin_run_script.m
@@ -27,15 +27,17 @@ EXECNAME = 'gyacomo23_dp'; % double precision
% run lin_DTT_HM_rho98
% run lin_DIIID_LM_rho90
% run lin_DIIID_LM_rho95
+% run lin_DIIID_LM_rho95_HEL
% run lin_JET_rho97
% run lin_Entropy
% run lin_ITG
% run lin_KBM
% run lin_RHT
-rho = 0.95; TRIANG = 'NT'; READPROF = 0;
-prof_folder = ['parameters/profiles/DIIID_Austin_et_al_2019/',TRIANG,'/'];
+% run lin_Ivanov
+rho = 0.95; TRIANG = 'NT'; READPROF = 1;
+% prof_folder = ['parameters/profiles/DIIID_Austin_et_al_2019/',TRIANG,'/'];
% prof_folder = ['parameters/profiles/DIIID_Oak_Nelson/',TRIANG,'/'];
-% prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_high_density/',TRIANG,'/'];
+prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_high_density/',TRIANG,'/'];
run lin_DIIID_data
% run lin_STEP_EC_HD_psi71
% run lin_STEP_EC_HD_psi49
@@ -43,47 +45,31 @@ if 1
% Plot the profiles
plot_params_vs_rho(geom,prof_folder,rho,0.5,1.1,Lref,mref,Bref,READPROF);
end
+% SIMID = ['rho_scan_DIIID_AUSTIN_2019/3x2x192x96x32/rho',num2str(rho)];
%% Change parameters
-% EXBRATE = 0.0; % Background ExB shear flow
-% K_Ti = 5.3;
-% NU = 0.001;
-CO = 'LD';
-GKCO = 1;
-kymin= 0.3; LY = 2*pi/kymin;
+% GEOMETRY = 's-alpha';
+DELTA =0.0;
+K_Ni = 0; K_Ne = 0;
+% DELTA = 0.0;
+% DELTA = 0.2;
+S_DELTA = DELTA/2;
+LY = 2*pi/0.25;
+TMAX = 20;
NY = 2;
-NX = 4;
-NZ = 32;
-PMAX = 2;
-JMAX = PMAX/2;
-DT = 20e-4;
-EXBRATE = 0;
-% EPS = 0.25;
-% KAPPA = 1.0;
-S_DELTA = min(2.0,S_DELTA);
-% DELTA = -DELTA;
-% PT parameters
-% TAU = 5.45;
-% K_Ne = 0; %vs 2.79
-% K_Ni = 0;
-% K_Te = 9.6455;%vs 17.3
-% K_Ti = 3.3640;%vs 5.15
-% BETA = 7.9e-4;
-% NU = 0.1;
-MU_X = 0.0; MU_Y = 0.0;
-SIGMA_E = 0.04;
-TMAX = 1;
-% DTSAVE0D = 200*DT;
-DTSAVE3D = 0.1;
-%%-------------------------------------------------------------------------
+DT = 0.01;
+TAU = 1; NU = 0.05;
+% TAU = 1e-3; K_Ti = K_Ti/2/TAU; NU = 3*NU/8/TAU; ADIAB_E = 1; NA = 1;
+% MU_X = 1; MU_Y = 1;
%% RUN
setup
+system(['cp ../results/',SIMID,'/',PARAMS,'/fort_00.90 ../results/',SIMID,'/.']);
% system(['rm fort*.90']);
% Run linear simulation
if RUN
MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
% RUN =['time ',mpirun,' -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
- RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
- % RUN =['time ',mpirun,' -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0;'];
+ RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0;'];
+ % RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 2 2 1 0;'];
% RUN =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 2 2 2 0;'];
% RUN =['time ',mpirun,' -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
% RUN = ['./../../../bin/gyacomo23_sp 0;'];
@@ -102,15 +88,15 @@ data = {}; % Initialize data as an empty cell array
% load grids, inputs, and time traces
data = compile_results_low_mem(data,LOCALDIR,J0,J1);
-if 0 % 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');
options.NORMALIZED = 0;
options.TIME = data.Ts3D;
% Time window to measure the growth of kx/ky modes
-options.KY_TW = [0.5 1.0]*data.Ts3D(end);
-options.KX_TW = [0.5 1.0]*data.Ts3D(end);
+options.KY_TW = [0.25 1.0]*data.Ts3D(end);
+options.KX_TW = [0.25 1.0]*data.Ts3D(end);
options.NMA = 1; % Set NMA option to 1
options.NMODES = 999; % Set how much modes we study
options.iz = 'avg'; % Compressing z
@@ -118,8 +104,38 @@ options.ik = 1; %
options.GOK2 = 0; % plot gamma/k^2
options.fftz.flag = 0; % Set fftz.flag option to 0
options.FIELD = 'phi';
-options.SHOWFIG = 1;
-[fig, kykx, wkykx, ekykx] = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
+options.SHOWFIG = 1;
+[fig, wkykx, ekykx] = mode_growth_meter(data,options);
+% %%
+% kx = (1:data.grids.Nx/2)'*2*pi/data.fort_00.GRID.Lx;
+ky = (1:data.grids.Ny/2)'*2*pi/data.fort_00.GRID.Ly;
+gkxky = real(wkykx(2:end,1:data.grids.Nx/2))';
+gkxky(isnan(gkxky)) =0;
+gkxky(isinf(gkxky)) =0;
+figure; plot(ky,gkxky(1,:));
+% gkxky(gkxky<0) =0;
+% % gkxky = imgaussfilt(gkxky,1);
+% %
+% wkxky = imag(wkykx(2:end,1:data.grids.Nx/2))';
+% wkxky(isnan(wkxky)) =0;
+% wkxky(isinf(wkxky)) =0;
+% % wkxky(wkxky<0) =0;
+% % wkxky = imgaussfilt(wkxky,1.5);
+% %
+% figure;
+% subplot(121)
+% contourf(kx,ky,gkxky',10)
+% % clim(0.5*[0 1]);
+% % colormap(bluewhitered); colorbar;
+% xlim([0.025 1]);
+% xlabel('$k_x\rho_s$'); ylabel('$k_y\rho_s$')
+% subplot(122)
+% contourf(kx,ky,wkxky',10)
+% % clim(1*[0 1]);
+% % colormap(bluewhitered); colorbar
+% xlim([0.025 1]);
+% xlabel('$k_x\rho_s$'); ylabel('$k_y\rho_s$')
+% % save_figure(data,fig,'.png')
end
if (0 && NZ>4)
@@ -144,13 +160,13 @@ ikx = 2; iky = 1; t0 = 1; t1 = data.Ts3D(end);
[~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D));
plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3));
t_ = data.Ts3D(it0:it1);
-% TH = 1.635*EPS^1.5 + 0.5*EPS^2+0.35*EPS^2.5; theory = 1/(1+Q0^2*TH/EPS^2);
+TH = 1.635*EPS^1.5 + 0.5*EPS^2+0.36*EPS^2.5; theory = 1/(1+Q0^2*TH/EPS^2);
clr_ = lines(20);
figure
-plot(t_, plt(data.PHI)); hold on;
+plot(t_, -plt(data.PHI)); hold on;
plot(t_,0.5* exp(-t_*NU)+theory,'--k');
plot([t_(1) t_(end)],theory*[1 1],'-k');
-plot([t_(1) t_(end)],mean(plt(data.PHI))*[1 1],'-g');
+plot([t_(1) t_(end)],-mean(plt(data.PHI))*[1 1],'-g');
xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.grids.kx(ikx),data.grids.ky(iky)))
end
@@ -160,7 +176,7 @@ if 0
plot_metric(data);
end
-if 1
+if 0
%% Compiled plot for lin EM analysis
[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
if data.inputs.BETA > 0
@@ -185,16 +201,22 @@ if 1
options.FIELD = 'phi';
options.SHOWFIG = 0;
[~, chi, phib, psib, ~] = plot_ballooning(data,options);
- [~, kykx, wkykx, ekykx] = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
+ [fig, wkykx, ekykx] = mode_growth_meter(data,options);
+ kx = (1:data.grids.Nx/2)'*2*pi/data.fort_00.GRID.Lx;
+ ky = (1:data.grids.Ny/2)'*2*pi/data.fort_00.GRID.Ly;
[~,~,R,Z] = plot_metric(data,options);
figure;
subplot(3,2,1); plot(chi,real(phib),'-b'); hold on;
plot(chi,imag(phib),'-r'); xlabel('$\chi$'); ylabel('$\phi$')
subplot(3,2,3); plot(chi,real(psib),'-b'); hold on;
plot(chi,imag(psib),'-r'); xlabel('$\chi$'); ylabel('$\psi$')
- subplot(3,2,5); plot(squeeze(kykx(:,1)),squeeze(real(wkykx(:,1))),'o--'); hold on;
- plot(squeeze(kykx(:,1)),squeeze(imag(wkykx(:,1))),'o--');
- xlabel('$k_y\rho_s$'); ylabel('$\gamma,\omega$');
+ subplot(3,2,5); errorbar(ky,squeeze(real(wkykx(2:end,1))),...
+ squeeze(real(ekykx(2:end,1))),...
+ 'ok--','MarkerSize',8,'LineWidth',1.5); hold on;
+ errorbar(ky,squeeze(imag(wkykx(2:end,1))),...
+ squeeze(imag(ekykx(2:end,1))),...
+ '*k--','MarkerSize',8,'LineWidth',1.5);
+ xlabel('$k_y\rho_s$'); ylabel('$\gamma (o),\,\omega (*)$');
R = R*geom.R0; Z = Z*geom.R0;
subplot(1,2,2); plot(R,Z,'-k'); xlabel('$R$'); ylabel('$Z$');
axis equal
diff --git a/wk/lin_scan_script.m b/wk/lin_scan_script.m
index bf36d44d2d156a2c46e53834f47edaa055c7b349..3883f0e9b4fff3be942ca2dc74e8f14c0b89dbbc 100644
--- a/wk/lin_scan_script.m
+++ b/wk/lin_scan_script.m
@@ -28,24 +28,29 @@ EXECNAME = 'gyacomo23_dp'; % double precision
% run lin_Entropy
% run lin_ITG
% run lin_RHT
-rho = 0.95; TRIANG = 'NT'; READPROF = 1;
+rho = 0.95; TRIANG = 'PT'; READPROF = 1;
% prof_folder = ['parameters/profiles/DIIID_Austin_et_al_2019/',TRIANG,'/'];
% prof_folder = ['parameters/profiles/DIIID_Oak_Nelson/',TRIANG,'/'];
prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_high_density/',TRIANG,'/'];
run lin_DIIID_data
%% Change parameters
-NU = 1;
-TAU = 1;
+% NU = 1;
+% TAU = 1;
NY = 2;
EXBRATE = 0;
+S_DELTA = min(2.0,S_DELTA);
SIGMA_E = 0.023;
+NEXC = 0;
+LX = 120;
%% Scan parameters
SIMID = [SIMID,'_scan'];
P_a = [2 4];
% P_a = 2;
-ky_a = [0.01 0.02 0.05 0.1 0.2 0.5 1.0 2.0 5.0 10.0];
-CO = 'LD';
+ky_a = [0.01 0.02 0.05 0.1 0.2 0.5 1.0 2.0 5.0 10.0];
+% ky_a = 4.0;
+dt_a = logspace(-2,-3,numel(ky_a));
+CO = 'DG';
%% Scan loop
% arrays for the result
g_ky = zeros(numel(ky_a),numel(P_a));
@@ -58,10 +63,10 @@ for PMAX = P_a
i = 1;
for ky = ky_a
LY = 2*pi/ky;
- DT = 2e-4;%/(1+log(ky/0.05));%min(1e-2,1e-3/ky);
- TMAX = 20;%min(10,1.5/ky);
- DTSAVE0D = 0.1;
- DTSAVE3D = 0.01;
+ DT = dt_a(i);%1e-3;%/(1+log(ky/0.05));%min(1e-2,1e-3/ky);
+ TMAX = DT*10000;%2;%min(10,1.5/ky);
+ DTSAVE0D = 100*DT;
+ DTSAVE3D = 10*DT;
%% RUN
setup
% naming
@@ -94,21 +99,27 @@ for PMAX = P_a
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 ...
- [wkykx,ekykx] = compute_growth_rates(data_.PHI(:,:,:,it1:it2),data_.Ts3D(it1:it2));
+ options.NORMALIZED = 0;
+ options.TIME = data_.Ts3D;
+ % Time window to measure the growth of kx/ky modes
+ options.KY_TW = [0.7 1.0]*data_.Ts3D(end);
+ options.KX_TW = [0.7 1.0]*data_.Ts3D(end);
+ options.NMA = 1; % Set NMA option to 1
+ options.NMODES = 999; % Set how much modes we study
+ options.iz = 'avg'; % Compressing z
+ options.ik = 1; %
+ options.GOK2 = 0; % plot gamma/k^2
+ options.fftz.flag = 0; % Set fftz.flag option to 0
+ options.FIELD = 'phi';
+ options.SHOWFIG = 0;
+ [fig, wkykx, ekykx] = mode_growth_meter(data_,options);
+ % [wkykx,ekykx] = compute_growth_rates(data_.PHI(:,:,:,it1:it2),data_.Ts3D(it1:it2));
g_ky (i,j) = real(wkykx(2,1));
g_std(i,j) = real(ekykx(2,1));
w_ky (i,j) = imag(wkykx(2,1));
w_std(i,j) = imag(ekykx(2,1));
[gmax, ikmax] = max(g_ky(i,j));
-
+
msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
end
i = i + 1;
diff --git a/wk/load_metadata_scan.m b/wk/load_metadata_scan.m
index a9df5c95be4a524979dbae0448bf4f1f2461c3c2..89db4acbe3035dde9e6d5f3226eba85105138ded 100644
--- a/wk/load_metadata_scan.m
+++ b/wk/load_metadata_scan.m
@@ -11,7 +11,7 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
% datafname = 'lin_DIIID_Oak_Nelson_high_density_NT_scan/6x32_ky_0.01_10_P_2_8_kN_1.0883_LDGK_0.0080915_be_0.0015991.mat';
% rho = 0.95
% datagname = 'lin_DIIID_Oak_Nelson_high_density_PT_scan/6x32_ky_0.01_10_P_2_8_kN_0.62888_LDGK_0.0046858_be_0.0048708.mat';
-datafname = 'lin_DIIID_Oak_Nelson_high_density_NT_scan/6x32_ky_0.01_10_P_2_8_kN_1.6989_LDGK_0.0167_be_0.00075879.mat';
+datafname = 'lin_DIIID_Oak_Nelson_high_density_PT_scan/6x32_ky_0.01_10_P_2_4_kN_0.62888_DGGK_0.0046858_be_0.0048708.mat';
%% Chose if we filter gamma>0.05
FILTERGAMMA = 1;
diff --git a/wk/old scripts/p3_geom_scan_analysis.m b/wk/old scripts/p3_geom_scan_analysis.m
new file mode 100644
index 0000000000000000000000000000000000000000..1c63315622f6a2a84a842307f9f66e86d04c9a7f
--- /dev/null
+++ b/wk/old scripts/p3_geom_scan_analysis.m
@@ -0,0 +1,306 @@
+% Get the current directory (wk)
+curdir = pwd;
+NCONTOUR = 0;
+% give ref value and param names
+REFVAL= 0;
+% normalize to the max all data
+NORM_ALL= 0;
+% normalize to the max each line
+NORM_LIN= 0;
+% normalize to the max each column
+NORM_COL= 0;
+% Get and plot the fluxsurface
+GETFLUXSURFACE = 0;
+
+% partition= '../results/paper_3/';
+% Get the scan directory
+switch 6
+ case 1 % delta kappa scan
+ casename = 'DTT rho85';
+ partition= '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/';
+ scandir = 'P2_J1_delta_kappa_scan'; scanname= '(2,1)';
+ % scandir = 'P4_J2_delta_kappa_scan'; scanname= '(4,2)';
+ nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0.23; scale1 =1.0;
+ nml2 = 'GEOMETRY'; pnam2 = '$\kappa$'; attr2 = 'kappa'; pref2 = 1.53; scale2 =1.0;
+ t1 = 50; t2 = 150;
+ case 2 % shear safety factor scan
+ casename = 'DTT rho85';
+ partition= '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/';
+ scandir = 'P2_J1_PT_sfact_shear_scan'; scanname= '(2,1)';
+ % scandir = 'P2_J1_NT_sfact_shear_scan'; scanname= '(2,1)';
+ nml1 = 'GEOMETRY'; pnam1 = '$\hat s$'; attr1 = 'shear'; pref1 = 3.63; scale1 =1.0;
+ nml2 = 'GEOMETRY'; pnam2 = '$q_0$'; attr2 = 'q0'; pref2 =-2.15; scale2 =1.0;
+ t1 = 50; t2 = 150;
+ case 3
+ casename = 'DTT rho85';
+ partition= '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/';
+ % scandir = 'P2_J1_delta_nuDGGK_scan'; scanname= 'DG (2,1)';
+ % scandir = 'P4_J2_delta_nuDGGK_scan'; scanname= 'DG (4,2)';
+ % scandir = 'P8_J4_delta_nuDGGK_conv_test'; scanname= 'DG (8,4)';
+ % scandir = 'P2_J1_delta_nuSGGK_scan'; scanname= 'SG (2,1)';
+ % scandir = 'P4_J2_delta_nuSGGK_scan'; scanname= 'SG (4,2)';
+ % scandir = 'P8_J4_delta_nuSGGK_conv_test'; scanname= 'SG (8,4)';
+ % scandir = 'P4_J2_delta_nuSGGKii_scan'; scanname= 'SGii (4,2)';
+ scandir = 'P2_J1_delta_nuLDGK_scan'; scanname= 'LD (2,1)';
+ % scandir = 'P4_J2_delta_nuLDGK_scan'; scanname= 'LD (4,2)';
+ nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0.23; scale1 =1.0;
+ nml2 = 'MODEL'; pnam2 = '$\nu$'; attr2 = 'nu'; pref2 = 0.5; scale2 =1.0;
+ t1 = 50; t2 = 150;
+ case 4 % shear delta scan
+ casename = 'DIIID rho95';
+ partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/';
+ scandir = '3x2x192x48x32_RT_2500_delta_shear_scan'; scanname= 'CI DG RT 2500';
+ % scandir = '3x2x256x64x48_delta_shear_scan';
+ nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0;
+ nml2 = 'GEOMETRY'; pnam2 = '$\hat s$'; attr2 = 'shear'; pref2 = 0.8; scale2 =1.0;
+ t1 = 50; t2 = 150;
+ case 5 % delta K_T tau=1
+ casename = 'DIIID rho95 $\tau=1$';
+ partition= '/misc/gyacomo23_outputs/paper_3/DIIID_tau_1_rho95_geom_scan/';
+ scandir = '3x2x192x48x32_delta_RT_scan'; scanname= '(2,1)';
+ % scandir = '5x3x192x48x32_delta_RT_scan'; scanname= '(4,2)';
+ % scandir = 'delta_RT_scan_PJ_21'; scanname= '(2,1)';
+ nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0;
+ nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 0.8; scale2 =1.0;
+ t1 = 200; t2 = 500;
+ case 6 % delta K_T cold ions
+ casename = 'DIIID rho95 $\tau=10^{-3}$';
+ partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/';
+ scandir = '3x2x192x48x32_delta_RT_scan'; scanname= '(2,1)';
+ nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0;
+ nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 0.8; scale2 =1e3/2;
+ t1 = 200; t2 = 480;
+ case 7 % delta s_delta
+ casename = 'DIIID rho95 $\tau=10^{-3}$';
+ partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/';
+ scandir = '3x2x192x48x32_RT_1000_delta_sdelta_scan'; scanname= 'RT=1000 (2,1)';
+ % scandir = '';
+ nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0;
+ nml2 = 'GEOMETRY'; pnam2 = '$s_\delta$'; attr2 = 's_delta'; pref2 = 0.8; scale2 =1.0;
+ t1 = 200; t2 = 295;
+ case 8 % eps q0
+ casename = 'DIIID rho95 $\tau=10^{-3}$';
+ partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/';
+ scandir = '3x2x192x48x32_RT_1000_eps_q0_scan/PT'; scanname= 'PT, RT=1000 (2,1)';
+ % scandir = '3x2x192x48x32_RT_1000_eps_q0_scan/NT'; scanname= 'NT, RT=1000 (2,1)';
+ % scandir = '';
+ nml1 = 'GEOMETRY'; pnam1 = '$\epsilon$'; attr1 = 'eps'; pref1 = 0; scale1 =1.0;
+ nml2 = 'GEOMETRY'; pnam2 = '$q_0$'; attr2 = 'q0'; pref2 = 0; scale2 =1.0;
+ t1 = 200; t2 = 400;
+ case 9 % CBC Dimits shift
+ casename = 'HEL CBC';
+ partition= '/misc/gyacomo23_outputs/paper_3/HEL_CBC/';
+ scandir = '128x32x24'; scanname= 'CBC HEL';
+ nml1 = 'SPECIES'; pnam1 = '$R_T$'; attr1 = 'k_T_'; pref1 = 0; scale1 =1.0;
+ nml2 = 'SPECIES'; pnam2 = '$R_N$'; attr2 = 'k_N_'; pref2 = 0; scale2 =1.0;
+ t1 = 1000; t2 = 2000;
+end
+scanname= [casename scanname];
+scandir = [partition,scandir,'/'];
+% Get a list of all items in the current directory
+contents = dir(scandir);
+
+% Iterate through the contents
+Qxavg = []; Qxerr = []; para1 = []; para2 = []; R = []; Z = [];
+Qxt = struct();
+for i = 1:length(contents)
+ % Check if the item is a directory and not '.' or '..'
+ if contents(i).isdir && ~strcmp(contents(i).name, '.') ...
+ && ~strcmp(contents(i).name, '..')
+ % Get and display the name of the subdirectory
+ subdir = [scandir,contents(i).name];
+ disp(['Subdirectory: ' contents(i).name]);
+ % Get parameters
+ param = read_namelist([subdir,'/fort_00.90']);
+ para1 = [para1 param.(nml1).(attr1)];
+ para2 = [para2 param.(nml2).(attr2)];
+ % Now you are in the subdirectory. You can perform operations here.
+ [t_all, Pxi_all, Qxi_all, Pxe_all, Qxe_all] = read_flux_out_XX(subdir);
+ if(numel(Qxe_all) > 1)
+ Qxtot = Qxi_all+Qxe_all;
+ else
+ Qxtot = Qxi_all;
+ end
+ Qxt.(['dat_',num2str(i)]) = struct();
+ Qxt.(['dat_',num2str(i)]).Qx = Qxtot;
+ Qxt.(['dat_',num2str(i)]).t = t_all;
+ Qxt.(['dat_',num2str(i)]).name = contents(i).name;
+ if(numel(t_all) > 1)
+ disp(num2str(t_all(end)))
+ [~,it1] = min(abs(t_all-t1));
+ [~,it2] = min(abs(t_all-t2));
+ steady_slice = it1:it2;
+ if(t_all(end) >= t2)
+ [fullAvg,sliceAvg,sliceErr] = sliceAverage(Qxtot(steady_slice),3);
+ Qxavg = [Qxavg fullAvg];
+ Qxerr = [Qxerr mean(sliceErr)];
+ else
+ Qxavg = [Qxavg nan];
+ Qxerr = [Qxerr nan];
+ end
+ else
+ Qxavg = [Qxavg nan];
+ Qxerr = [Qxerr nan];
+ end
+ end
+ if GETFLUXSURFACE
+ data = load([subdir,'/RZ.txt']);
+ R_ = data(:, 1);
+ Z_ = data(:, 2);
+ R_ = [R_;R_(1)]'; Z_ = [Z_;Z_(1)]';
+ R = [R ; R_]; Z = [Z ; Z_];
+ end
+
+end
+if 0
+%% plot time traces
+attr = fieldnames(Qxt);
+Nsim = numel(attr);
+figure
+% compute growth at the begining
+tw = [10 40];
+gr = 1:Nsim; err = 1:Nsim;
+for i = 1:1:Nsim
+ tmp_ = Qxt.(attr{i});
+ t = tmp_.t;
+ y = tmp_.Qx;
+ plot(t,y,'DisplayName',tmp_.name); hold on;
+ [~,it1] = min(abs(t-tw(1)));
+ [~,it2] = min(abs(t-tw(2)));
+ [gr_, err_] = compute_growth(t(it1:it2),y(it1:it2));
+ gr(i) = gr_; err(i) = err_;
+end
+%%
+toplot = real(reshape(gr,sz))';
+toplot = toplot(idx1,idx2);
+
+figure
+imagesc_custom(xx_,yy_,toplot); hold on
+end
+%% reshaping, sorting and plotting
+p1 = unique(para1)/scale1;
+p2 = unique(para2)/scale2;
+N1 = numel(p1);
+N2 = numel(p2);
+
+if para1(1) == para1(2)
+ sz = [N2 N1];
+ TRANSPOSE = 1;
+else
+ sz = [N1 N2];
+ TRANSPOSE = 0;
+end
+
+Zavg = reshape(Qxavg,sz);
+Zerr = reshape(Qxerr,sz);
+XX = reshape(para1/scale1,sz);
+YY = reshape(para2/scale2,sz);
+
+if TRANSPOSE
+ Zavg = Zavg';
+ Zerr = Zerr';
+ XX = XX';
+ YY = YY';
+end
+
+[~,idx1] = sort(XX(:,1));
+[~,idx2] = sort(YY(1,:));
+Zavg = Zavg(idx1,idx2);
+Zerr = Zerr(idx1,idx2);
+XX = XX(idx1,idx2);
+YY = YY(idx1,idx2);
+
+% compute the
+if REFVAL
+ if NORM_ALL
+ Qxname = '$\bar Q_{tot}/\bar Q_{max}[\%]$';
+ [tmp,iref1] = max(Zavg);
+ [~, iref2] = max(tmp);
+ iref1 = iref1(iref2);
+ else
+ Qxname = '$\langle (Q_{tot}-Q_{ref})/Q_{ref} \rangle_t[\%]$';
+ if pref1 ~= 999
+ [~,iref1] = min(abs(XX(:,1)-pref1));
+ else
+ iref1 = 1:N1;
+ end
+ if pref2 ~= 999
+ [~,iref2] = min(abs(YY(1,:)-pref2));
+ else
+ iref2 = 1:N2;
+ end
+ end
+ iref1 = ones(N1,1).*iref1;
+ iref2 = ones(N2,1).*iref2;
+ xref = XX(iref1,iref2);
+ yref = YY(iref1,iref2);
+ Qxref = Zavg(iref1,iref2);
+ Qrefname = ['$Q_{ref}=$',num2str(Qxref(1,1))];
+else
+ Qref = 1;
+ if NORM_LIN
+ Qxname = '$\bar Q_{tot}/\bar Q_{max}[\%]$, per line';
+ for il = 1:sz(1)
+ maxline = max(Zavg(:,il));
+ Zavg(:,il) = Zavg(:,il)./maxline;
+ Zerr(:,il) = Zerr(:,il)./maxline;
+ end
+ elseif NORM_COL
+ Qxname = '$\bar Q_{tot}/\bar Q_{max}[\%]$, per column';
+ for ic = 1:sz(2)
+ maxcol = max(Zavg(ic,:));
+ Zavg(ic,:) = Zavg(ic,:)./maxcol;
+ Zerr(ic,:) = Zerr(ic,:)./maxcol;
+ end
+ else
+ Qxname = '$\langle Q_{tot} \rangle_t$';
+ end
+end
+
+% Figure
+figure
+subplot(1,2,1)
+[xx_,yy_] = meshgrid(XX(:,1),YY(1,:));
+if REFVAL
+ if NORM_ALL || NORM_COL || NORM_LIN
+ toplot = (Zavg./Qxref)'
+ CLIM = [0 1];
+ else
+ toplot = ((Zavg-Qxref)./Qxref * 100)';
+ CLIM = 'auto';
+ end
+else
+ toplot = Zavg';
+ CLIM = 'auto';
+end
+if NCONTOUR <= 0
+ imagesc_custom(xx_,yy_,toplot); hold on
+else
+ contourf(XX(:,1),YY(1,:),Zavg',NCONTOUR); hold on
+end
+if REFVAL && ~((pref1==999) || (pref2==999))
+ plot(xref(1,1),yref(1,1),'xk','MarkerSize',14,'DisplayName',Qrefname)
+ legend('show')
+end
+xlabel(pnam1); ylabel(pnam2);
+title(scanname)
+colormap(bluewhitered); colorbar; clim(CLIM);
+if ~REFVAL
+ colormap(jet);
+end
+subplot(1,2,2)
+clrs = jet(N2);
+for i = 1:N2
+ errorbar(XX(:,i),Zavg(:,i),Zerr(:,i),...
+ 'DisplayName',[pnam2,'=',num2str(p2(i))],...
+ 'Color',clrs(i,:));
+ hold on;
+end
+if REFVAL && ~((pref1==999) || (pref2==999))
+ plot(xref(1,1),0,'xk','MarkerSize',14,'DisplayName',Qrefname)
+end
+grid on
+xlabel(pnam1); ylabel('$\langle Q_{tot} \rangle_t$');
+legend('show','Location','northwest');
+title([param.COLLISION.collision_model{1}, ...
+ ', $(P,J)=(',num2str(param.GRID.pmax),',',num2str(param.GRID.jmax),')$'])
diff --git a/wk/p3_geom_scan_analysis.m b/wk/p3_geom_scan_analysis.m
index f364ab5616f557504e4d7fb015fb9a3fedd702bf..19111f49919affb0cdcee5ec77c672add848670c 100644
--- a/wk/p3_geom_scan_analysis.m
+++ b/wk/p3_geom_scan_analysis.m
@@ -14,86 +14,37 @@ GETFLUXSURFACE = 0;
% partition= '../results/paper_3/';
% Get the scan directory
-switch 5
- case 1 % delta kappa scan
- casename = 'DTT rho85';
- partition= '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/';
- scandir = 'P2_J1_delta_kappa_scan'; scanname= '(2,1)';
- % scandir = 'P4_J2_delta_kappa_scan'; scanname= '(4,2)';
- nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0.23; scale1 =1.0;
- nml2 = 'GEOMETRY'; pnam2 = '$\kappa$'; attr2 = 'kappa'; pref2 = 1.53; scale2 =1.0;
- t1 = 50; t2 = 150;
- case 2 % shear safety factor scan
- casename = 'DTT rho85';
- partition= '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/';
- scandir = 'P2_J1_PT_sfact_shear_scan'; scanname= '(2,1)';
- % scandir = 'P2_J1_NT_sfact_shear_scan'; scanname= '(2,1)';
- nml1 = 'GEOMETRY'; pnam1 = '$\hat s$'; attr1 = 'shear'; pref1 = 3.63; scale1 =1.0;
- nml2 = 'GEOMETRY'; pnam2 = '$q_0$'; attr2 = 'q0'; pref2 =-2.15; scale2 =1.0;
- t1 = 50; t2 = 150;
- case 3
- casename = 'DTT rho85';
- partition= '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/';
- % scandir = 'P2_J1_delta_nuDGGK_scan'; scanname= 'DG (2,1)';
- % scandir = 'P4_J2_delta_nuDGGK_scan'; scanname= 'DG (4,2)';
- % scandir = 'P8_J4_delta_nuDGGK_conv_test'; scanname= 'DG (8,4)';
- % scandir = 'P2_J1_delta_nuSGGK_scan'; scanname= 'SG (2,1)';
- % scandir = 'P4_J2_delta_nuSGGK_scan'; scanname= 'SG (4,2)';
- % scandir = 'P8_J4_delta_nuSGGK_conv_test'; scanname= 'SG (8,4)';
- % scandir = 'P4_J2_delta_nuSGGKii_scan'; scanname= 'SGii (4,2)';
- scandir = 'P2_J1_delta_nuLDGK_scan'; scanname= 'LD (2,1)';
- % scandir = 'P4_J2_delta_nuLDGK_scan'; scanname= 'LD (4,2)';
- nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0.23; scale1 =1.0;
- nml2 = 'MODEL'; pnam2 = '$\nu$'; attr2 = 'nu'; pref2 = 0.5; scale2 =1.0;
- t1 = 50; t2 = 150;
- case 4 % shear delta scan
- casename = 'DIIID rho95';
- partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/';
- scandir = '3x2x192x48x32_RT_2500_delta_shear_scan'; scanname= 'CI DG RT 2500';
- % scandir = '3x2x256x64x48_delta_shear_scan';
- nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0;
- nml2 = 'GEOMETRY'; pnam2 = '$\hat s$'; attr2 = 'shear'; pref2 = 0.8; scale2 =1.0;
- t1 = 50; t2 = 150;
- case 5 % delta K_T tau=1
+switch 2
+ case 1 % delta K_T tau=1
casename = 'DIIID rho95 $\tau=1$';
partition= '/misc/gyacomo23_outputs/paper_3/DIIID_tau_1_rho95_geom_scan/';
- scandir = '3x2x192x48x32_delta_RT_scan'; scanname= '(2,1)';
- % scandir = '5x3x192x48x32_delta_RT_scan'; scanname= '(4,2)';
+ % scandir = '3x2x192x48x32_nu_0.05_delta_RT_scan'; scanname= '(2,1)';
+ % scandir = '3x2x192x48x32_nu_0.1_delta_RT_scan'; scanname= '(2,1)';
+ % scandir = '3x2x192x48x24_nu_0.1_delta_RT_scan'; scanname= '(2,1)';
+ % scandir = '3x2x192x48x32_nu_1.0_delta_RT_scan'; scanname= '(2,1)';
+ % scandir = '5x3x192x48x32_nu_0.05_delta_RT_scan'; scanname= '(4,2)';
+ scandir = '5x3x192x48x32_nu_1.0_delta_RT_scan'; scanname= '(4,2)';
+ % scandir = '5x3x192x48x32_delta_RT_scan'; scanname= '(2,1)';
% scandir = 'delta_RT_scan_PJ_21'; scanname= '(2,1)';
nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0;
- nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 0.8; scale2 =1.0;
- t1 = 200; t2 = 800;
- case 6 % delta K_T cold ions
+ nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 5; scale2 =1.0;
+ t1 = 300; t2 = 500; zfactor = 1;
+ case 2 % delta K_T cold ions
casename = 'DIIID rho95 $\tau=10^{-3}$';
partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/';
- scandir = '3x2x192x48x32_delta_RT_scan'; scanname= '(2,1)';
- nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0;
- nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 0.8; scale2 =1e3;
- t1 = 200; t2 = 480;
- case 7 % delta s_delta
- casename = 'DIIID rho95 $\tau=10^{-3}$';
- partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/';
- scandir = '3x2x192x48x32_RT_1000_delta_sdelta_scan'; scanname= 'RT=1000 (2,1)';
- % scandir = '';
+ % scandir = '3x2x192x48x32_nu_0_delta_RT_scan'; scanname= '(2,1)';
+ scandir = '3x2x192x48x32_nu_0.05_delta_RT_scan'; scanname= '(2,1)';
+ nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1;
+ nml2 = 'SPECIES'; pnam2 = '$\kappa_T$'; attr2 = 'K_T_'; pref2 = 5; scale2 =500;
+ t1 = 80; t2 = 400; zfactor = 2;
+ case 3 % delta K_T HEL, better resolution
+ casename = 'DIIID rho95 $\tau=1$';
+ % partition= '/misc/gyacomo23_outputs/paper_3/geom_scan_DIIID_HEL/NU_50/';
+ partition= '/misc/gyacomo23_outputs/paper_3/geom_scan_DIIID_HEL/NU_20/';
+ scandir = '.'; scanname= 'CBC HEL';
nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0;
- nml2 = 'GEOMETRY'; pnam2 = '$s_\delta$'; attr2 = 's_delta'; pref2 = 0.8; scale2 =1.0;
- t1 = 200; t2 = 295;
- case 8 % eps q0
- casename = 'DIIID rho95 $\tau=10^{-3}$';
- partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/';
- scandir = '3x2x192x48x32_RT_1000_eps_q0_scan/PT'; scanname= 'PT, RT=1000 (2,1)';
- % scandir = '3x2x192x48x32_RT_1000_eps_q0_scan/NT'; scanname= 'NT, RT=1000 (2,1)';
- % scandir = '';
- nml1 = 'GEOMETRY'; pnam1 = '$\epsilon$'; attr1 = 'eps'; pref1 = 0; scale1 =1.0;
- nml2 = 'GEOMETRY'; pnam2 = '$q_0$'; attr2 = 'q0'; pref2 = 0; scale2 =1.0;
- t1 = 200; t2 = 400;
- case 9 % CBC Dimits shift
- casename = 'HEL CBC';
- partition= '/misc/gyacomo23_outputs/paper_3/HEL_CBC/';
- scandir = '128x32x24'; scanname= 'CBC HEL';
- nml1 = 'SPECIES'; pnam1 = '$R_T$'; attr1 = 'k_T_'; pref1 = 0; scale1 =1.0;
- nml2 = 'SPECIES'; pnam2 = '$R_N$'; attr2 = 'k_N_'; pref2 = 0; scale2 =1.0;
- t1 = 1000; t2 = 2000;
+ nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 5; scale2 =500;
+ t1 = 100; t2 = 150; zfactor = 1;
end
scanname= [casename scanname];
scandir = [partition,scandir,'/'];
@@ -115,11 +66,16 @@ for i = 1:length(contents)
para1 = [para1 param.(nml1).(attr1)];
para2 = [para2 param.(nml2).(attr2)];
% Now you are in the subdirectory. You can perform operations here.
- [t_all, Pxi_all, Qxi_all, Pxe_all, Qxe_all] = read_flux_out_XX(subdir);
+ out = read_flux_out_XX(subdir);
+ t_all = out.t;
+ Pxi_all = out.Pxi;
+ Qxi_all = out.Qxi;
+ Pxe_all = out.Pxe;
+ Qxe_all = out.Qxe;
if(numel(Qxe_all) > 1)
- Qxtot = Qxi_all+Qxe_all;
+ Qxtot = zfactor*(Qxi_all+Qxe_all);
else
- Qxtot = Qxi_all;
+ Qxtot = zfactor*(Qxi_all);
end
Qxt.(['dat_',num2str(i)]) = struct();
Qxt.(['dat_',num2str(i)]).Qx = Qxtot;
@@ -139,7 +95,7 @@ for i = 1:length(contents)
Qxerr = [Qxerr nan];
end
else
- Qxavg = [Qxavg nan];
+ Qxavg = [Qxavg nan];
Qxerr = [Qxerr nan];
end
end
@@ -158,7 +114,7 @@ attr = fieldnames(Qxt);
Nsim = numel(attr);
figure
% compute growth at the begining
-tw = [10 40];
+tw = [5 20];
gr = 1:Nsim; err = 1:Nsim;
for i = 1:1:Nsim
tmp_ = Qxt.(attr{i});
@@ -276,7 +232,7 @@ end
if NCONTOUR <= 0
imagesc_custom(xx_,yy_,toplot); hold on
else
- contourf(XX(:,1),YY(1,:),Zavg',NCONTOUR); hold on
+ contour(XX(:,1),YY(1,:),Zavg'); hold on
end
if REFVAL && ~((pref1==999) || (pref2==999))
plot(xref(1,1),yref(1,1),'xk','MarkerSize',14,'DisplayName',Qrefname)
@@ -304,3 +260,40 @@ xlabel(pnam1); ylabel('$\langle Q_{tot} \rangle_t$');
legend('show','Location','northwest');
title([param.COLLISION.collision_model{1}, ...
', $(P,J)=(',num2str(param.GRID.pmax),',',num2str(param.GRID.jmax),')$'])
+
+if 0
+%% plot minimum
+idxmax = 1:numel(Zavg(1,:));
+idxmin = 1:numel(Zavg(1,:));
+xmax = 1:numel(Zavg(1,:));
+xmin = 1:numel(Zavg(1,:));
+ymax = 1:numel(Zavg(1,:));
+ymin = 1:numel(Zavg(1,:));
+err = 1:numel(Zavg(1,:));
+
+x = linspace(min(p1),max(p1),128);
+for i=1:numel(Zavg(1,:))
+ [fit, dat] = polyfit(p1,Zavg(:,i)+0*Zerr(:,i),2);
+ [ymax(i),idx] = min(polyval(fit,x));
+ xmax(i) = x(idx);
+ [fit, dat] = polyfit(p1,Zavg(:,i)-0*Zerr(:,i),2);
+ [ymin(i),idx] = min(polyval(fit,x));
+ xmin(i) = x(idx);
+
+ [zmin,idx] = min(Zavg(:,i));
+ % err(i) = abs(zmin-Zavg(idx+1,i))/abs(zmin)+abs(zmin-Zavg(idx-1,i))/abs(zmin);
+end
+err = min(err,1);
+xavg = 0.5*(xmax+xmin);
+xerr = 0.5*abs(xmax-xmin);
+
+fit = polyfit(p2,xavg,1);
+y = linspace(min(p2),max(p2),128);
+
+figure
+plot(xavg+xerr,p2); hold on
+plot(xavg-xerr,p2); hold on
+plot(polyval(fit,y),y)
+plot(polyval(fit,y),y)
+plot(polyval(fit,y),y)
+end
\ No newline at end of file
diff --git a/wk/parallel_scaling.m b/wk/parallel_scaling.m
index ed35a92c497eedd0683c4d59d9a1d50616a051cb..5df3ebe177b84e9856bb87daadc5ea557ab760d8 100644
--- a/wk/parallel_scaling.m
+++ b/wk/parallel_scaling.m
@@ -1,11 +1,12 @@
-DIR = '/misc/gyacomo23_outputs/scaling/strong/17x9x256x256x32/'; scaling='strong';
+% DIR = '/misc/gyacomo23_outputs/scaling/strong/17x9x256x256x32/'; scaling='strong';
% DIR = '/misc/gyacomo23_outputs/scaling/strong/9x5x768x192x32/'; scaling='strong';
-% DIR = '/misc/gyacomo23_outputs/scaling/strong/3x2x768x192x32/'; scaling='strong';
+DIR = '/misc/gyacomo23_outputs/scaling/strong/3x2x768x192x32/'; scaling='strong';
% DIR = '/misc/gyacomo23_outputs/scaling/weak/Np_5x2x128x64x32/'; scaling='weak';
% DIR = '/misc/gyacomo23_outputs/scaling/weak/Ny_3x2x768x8x32/'; scaling='weak';
% DIR = '/misc/gyacomo23_outputs/scaling/weak/Nz_3x2x768x32x8/'; scaling='weak';
% DIR = '/misc/gyacomo23_outputs/scaling/weak/Nz_5x2x128x64x8/'; scaling='weak';
-% DIR = '/misc/gyacomo23_outputs/scaling/weak/Npyz_4x2x32x16x16/'; scaling='weak';
+% DIR = '/misc/gyacomo23_outputs/scaling/weak/Npyz_4x2x32x16x8/'; scaling='weak';
+% DIR = '/misc/gyacomo23_outputs/scaling/weak/Npyz_8x2x32x32x4/'; scaling='weak';
% Get a list of all items in the current directory
contents = dir(DIR);
@@ -76,31 +77,52 @@ Rt_ref = Rt_avg(Np_tot==1);
Np_ref = 1;
MaxN = max(Np_tot);
Np1Nz1 = logical((Np==1).*(Nz==1));
-Np2Nz1 = logical((Np==2).*(Nz==1));
Np1Nz2 = logical((Np==1).*(Nz==2));
Np1Nz4 = logical((Np==1).*(Nz==4));
+Np2Nz1 = logical((Np==2).*(Nz==1));
+Np2Nz2 = logical((Np==2).*(Nz==2));
+Np2Nz4 = logical((Np==2).*(Nz==4));
+Np4Nz1 = logical((Np==4).*(Nz==1));
+Np4Nz2 = logical((Np==4).*(Nz==2));
+Np4Nz4 = logical((Np==4).*(Nz==4));
%%
-figure; hold on;
+figure; hold on; msz = 8;
xlabel('Number of cores');
+clr_ = lines(3);
switch scaling
case 'strong'
- plot(Np_tot(Np1Nz1),Rt_ref./Rt_avg(Np1Nz1),...
- 'o','DisplayName','$N_{pp}=1, N_{pz}=1$');
- plot(Np_tot(Np2Nz1),Rt_ref./Rt_avg(Np2Nz1), ...
- 'o','DisplayName','$N_{pp}=2, N_{pz}=1$');
- plot(Np_tot(Np1Nz2),Rt_ref./Rt_avg(Np1Nz2), ...
- 'o','DisplayName','$N_{pp}=1, N_{pz}=2$');
- plot(Np_tot(Np1Nz4),Rt_ref./Rt_avg(Np1Nz4), ...
- 'o','DisplayName','$N_{pp}=1, N_{pz}=4$');
- plot(1:MaxN,(1:MaxN)/Np_ref,'--k','DisplayName','Ideal');
+ plot(Np_tot(Np1Nz1),Rt_ref./Rt_avg(Np1Nz1),'o', ...
+ 'Color',clr_(1,:),'DisplayName','$N_{pp}=1, N_{pz}=1$','MarkerSize',msz);
+ plot(Np_tot(Np1Nz2),Rt_ref./Rt_avg(Np1Nz2),'o',...
+ 'Color',clr_(2,:),'DisplayName','$N_{pp}=1, N_{pz}=2$','MarkerSize',msz);
+ plot(Np_tot(Np1Nz4),Rt_ref./Rt_avg(Np1Nz4),'o', ...
+ 'Color',clr_(3,:),'DisplayName','$N_{pp}=1, N_{pz}=4$','MarkerSize',msz);
+ plot(Np_tot(Np2Nz1),Rt_ref./Rt_avg(Np2Nz1),'x', ...
+ 'Color',clr_(1,:),'DisplayName','$N_{pp}=2, N_{pz}=1$','MarkerSize',msz);
+ plot(Np_tot(Np2Nz2),Rt_ref./Rt_avg(Np2Nz2),'x', ...
+ 'Color',clr_(2,:),'DisplayName','$N_{pp}=2, N_{pz}=2$','MarkerSize',msz);
+ plot(Np_tot(Np2Nz4),Rt_ref./Rt_avg(Np2Nz4),'x', ...
+ 'Color',clr_(3,:),'DisplayName','$N_{pp}=2, N_{pz}=4$','MarkerSize',msz);
+ plot(Np_tot(Np4Nz1),Rt_ref./Rt_avg(Np4Nz1),'+', ...
+ 'Color',clr_(1,:),'DisplayName','$N_{pp}=4, N_{pz}=1$','MarkerSize',msz);
+ plot(Np_tot(Np4Nz2),Rt_ref./Rt_avg(Np4Nz2),'+', ...
+ 'Color',clr_(2,:),'DisplayName','$N_{pp}=4, N_{pz}=2$','MarkerSize',msz);
+ plot(Np_tot(Np4Nz4),Rt_ref./Rt_avg(Np4Nz4),'+', ...
+ 'Color',clr_(3,:),'DisplayName','$N_{pp}=4, N_{pz}=4$','MarkerSize',msz);
+ plot(Np_tot,Rt_ref./Rt_avg,...
+ '.k','DisplayName','Strong scaling');
+ plot([1 1000],[1 1000],'--k','DisplayName','Ideal');
+ axis equal
set(gca,'XScale','log'); set(gca,'YScale','log');
ylabel('Speedup');
% title('Strong scaling speedup')
case 'weak'
hold on;
- filt = logical((Np==2).*(Ny>=1).*(Nz>=1));
- plot(Np_tot(filt),Rt_ref./Rt_avg(filt),...
- 'o','DisplayName','Effective');
+ % for n = [1 2 4 8 12]
+ filt = logical((Np>=1).*(Ny>=1).*(Nz>=1));
+ plot(Np_tot(filt),Rt_ref./Rt_avg(filt),...
+ 'o','DisplayName','Effective');
+ % end
plot(1:MaxN,ones(1,MaxN),'--k','DisplayName','Ideal');
ylim([0 1]);
set(gca,'XScale','log')
diff --git a/wk/parameters/lin_DIIID_LM_rho95.m b/wk/parameters/lin_DIIID_LM_rho95.m
index 2bd5ba9b2f383bfb8a885ecb30666f6b79194743..17c6cbf8552518eb05c6cbd2b62d9671e59d0730 100644
--- a/wk/parameters/lin_DIIID_LM_rho95.m
+++ b/wk/parameters/lin_DIIID_LM_rho95.m
@@ -96,4 +96,5 @@ BCKGD0 = 0.0e-5; % Initial background
k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; % Magnetic curvature strength
COLL_KCUT = 1; % Cutoff for collision operator
-ADIAB_I = 0; % adiabatic ion model
\ No newline at end of file
+ADIAB_I = 0; % adiabatic ion model
+EXBRATE = 0;
\ No newline at end of file
diff --git a/wk/parameters/lin_DIIID_LM_rho95_HEL.m b/wk/parameters/lin_DIIID_LM_rho95_HEL.m
new file mode 100644
index 0000000000000000000000000000000000000000..cbd940901c56c65a2f11ccf80320324342097e03
--- /dev/null
+++ b/wk/parameters/lin_DIIID_LM_rho95_HEL.m
@@ -0,0 +1,100 @@
+%% Reference values
+% See Neiser et al. 2019 Gyrokinetic GENE simulations of DIII-D near-edge L-mode plasmas
+%% Set simulation parameters
+SIMID = 'lin_DIIID_HM_rho90'; % Name of the simulation
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NU = 50; %(0.00235 in GENE)
+TAU = 0.001; % i/e temperature ratio
+K_Ne = 0; % ele Density '''
+K_Te = 0; % ele Temperature '''
+K_Ni = 0; % ion Density gradient drive
+K_Ti = 5.2256/2/TAU; % ion Temperature '''
+SIGMA_E = 0.0233380/sqrt(2); % 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; % electron plasma beta in prct
+MHD_PD = 1;
+%% Set up grid parameters
+P = 2;
+J = P/2;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 8; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
+LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.3; % Size of the squared frequency domain in y direction
+NZ = 32; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+
+%% GEOMETRY
+% GEOMETRY= 's-alpha';
+GEOMETRY= 'miller';
+Q0 = 4.8; % safety factor
+SHEAR = 2.55; % magnetic shear
+EPS = 0.3; % inverse aspect ratio
+KAPPA = 1.57; % elongation
+S_KAPPA = 0.48;
+DELTA = 0.2; % triangularity
+S_DELTA = 0.1;
+ZETA = 0.0; % squareness
+S_ZETA = 0.0;
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+PB_PHASE = 0;
+%% TIME PARAMETERS
+TMAX = 15; % Maximal time unit
+DT = 1e-2; % Time step
+DTSAVE0D = 0.5; % Sampling time for 0D arrays
+DTSAVE2D = -1; % Sampling time for 2D arrays
+DTSAVE3D = 0.5; % Sampling time for 3D arrays
+DTSAVE5D = 100; % Sampling time for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 1; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
+
+%% OUTPUTS
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.0; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+MU_Z = 5.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 1.0e-5; % Initial noise amplitude
+BCKGD0 = 0.0e-5; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1; % Cutoff for collision operator
+ADIAB_I = 0; % adiabatic ion model
+EXBRATE = 0;
\ No newline at end of file
diff --git a/wk/parameters/lin_ITG.m b/wk/parameters/lin_ITG.m
index 07d85ebc3aa01595ebd7e924280e8eaf25d875c9..d1c06b664d3094daa98e7e1c8e4c327e342cf16f 100644
--- a/wk/parameters/lin_ITG.m
+++ b/wk/parameters/lin_ITG.m
@@ -15,10 +15,8 @@ ADIAB_E = (NA==1); % adiabatic electron model
BETA = 0.0; % electron plasma beta
EXBRATE = 0.0; % Background ExB shear flow
%% Set up grid parameters
-P = 4;
-J = 2;%P/2;
-PMAX = P; % Hermite basis size
-JMAX = J; % Laguerre basis size
+PMAX = 4; % Hermite basis size
+JMAX = 2; % Laguerre basis size
NX = 8; % real space x-gridpoints
NY = 12; % real space y-gridpoints
LX = 2*pi/0.05; % Size of the squared frequency domain in x direction
diff --git a/wk/parameters/lin_Ivanov.m b/wk/parameters/lin_Ivanov.m
index a638d755eca8e64b372504c7a5f5d6abb581518c..c897d709f35e8a091b029f1ba0165886907193a3 100644
--- a/wk/parameters/lin_Ivanov.m
+++ b/wk/parameters/lin_Ivanov.m
@@ -4,11 +4,11 @@ SIMID = 'lin_Ivanov'; % Name of the simulation
%% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
TAU = 1e-3; % e/i temperature ratio
-NU = 0.1/TAU; % Collision frequency
+NU = 0.1*3/8/TAU; % Collision frequency
K_Ne = 0*2.22; % ele Density
K_Te = 0*6.96; % ele Temperature
K_Ni = 0*2.22; % ion Density gradient drive
-K_Ti = 0.36*2/TAU; % ion Temperature
+K_Ti = 1*2/TAU; % 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
@@ -51,8 +51,8 @@ JOB2LOAD = -1; % Start a new simulation serie
%% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
-CO = 'SG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-GKCO = 0; % Gyrokinetic operator
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 1; % Gyrokinetic operator
ABCO = 1; % INTERSPECIES collisions
INIT_ZF = 0; % Initialize zero-field quantities
% HRCY_CLOS = 'max_degree'; % Closure model for higher order moments
@@ -93,3 +93,10 @@ BCKGD0 = 0.0; % Initial background
k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; % Magnetic curvature strength
COLL_KCUT = 1.0; % Cutoff for collision operator
+S_KAPPA = 0.0;
+S_DELTA = 0.0;
+S_ZETA = 0.0;
+PB_PHASE= 0;
+ADIAB_I = 0;
+MHD_PD = 0;
+EXBRATE = 0;
\ No newline at end of file
diff --git a/wk/parameters/lin_RHT.m b/wk/parameters/lin_RHT.m
index 72137f75d6cb2f219d7ce3bf3bd81ecc77aea40e..0f39498d342d28c11b63c43d64e70132d1f07e0a 100644
--- a/wk/parameters/lin_RHT.m
+++ b/wk/parameters/lin_RHT.m
@@ -15,15 +15,15 @@ ADIAB_E = (NA==1); % adiabatic electron model
BETA = 0.0; % electron plasma beta
EXBRATE = 0.0; % Background ExB shear flow
%% Set up grid parameters
-P = 4;
-J = 2;%P/2;
+P = 256;
+J = 16;%P/2;
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size
NX = 2; % real space x-gridpoints
NY = 2; % real space y-gridpoints
LX = 2*pi/0.05; % Size of the squared frequency domain in x direction
LY = 2*pi/0.01; % Size of the squared frequency domain in y direction
-NZ = 8; % number of perpendicular planes (parallel grid)
+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))
@@ -42,7 +42,7 @@ NPOL = 1; % Number of poloidal turns
%% TIME PARAMETERS
TMAX = 50; % Maximal time unit
-DT = 1e-3; % Time step
+DT = 5e-3; % Time step
DTSAVE0D = 1; % Sampling time for 0D arrays
DTSAVE2D = -1; % Sampling time for 2D arrays
DTSAVE3D = 0.1; % Sampling time for 3D arrays
diff --git a/wk/particle_trajectory.m b/wk/particle_trajectory.m
new file mode 100644
index 0000000000000000000000000000000000000000..62da611a8b43f50d5b00184918e8d22416b0b2ed
--- /dev/null
+++ b/wk/particle_trajectory.m
@@ -0,0 +1,36 @@
+R = 1;
+a = 0.05;
+d = 0.02;
+Wc= 80;
+
+
+figure
+
+% Magnetic field line
+Bx = @(R,th) R*cos(th);
+By = @(R,th) R*sin(th);
+Bz = @(R,th) 0*(th);
+
+thB = linspace(0,2*pi,128);
+plot3(Bx(R,thB),By(R,thB),Bz(R,thB)); hold on;
+
+% Particle trochoid
+Tx = @(R,th) a*cos(th).*cos(Wc*th);
+Ty = @(R,th) a*sin(th).*cos(Wc*th);
+Tz = @(R,th) a*sin(Wc*th)+d*th;
+
+thP = linspace(0,2.5*pi,2048);
+Xp = Bx(R,thP)+Tx(R,thP);
+Yp = By(R,thP)+Ty(R,thP);
+Zp = Bz(R,thP)+Tz(R,thP);
+plot3(Xp,Yp,Zp)
+plot3(Xp(1),Yp(1),Zp(1),'xk','MarkerSize',10)
+plot3(Xp(end),Yp(end),Zp(end),'ok','MarkerSize',8,'MarkerFaceColor','k')
+
+
+% finish the plot
+axis equal
+xlim(1.2*R*[-1 1])
+ylim(1.2*R*[-1 1])
+zlim(0.25*[-1 1])
+axis off
\ No newline at end of file
diff --git a/wk/plot_params_vs_rho.m b/wk/plot_params_vs_rho.m
index 0b3d64860fbba428210049054fa658e9b577785f..0f12e391a80a4bc72a373a40bc098cdb4bec93ae 100644
--- a/wk/plot_params_vs_rho.m
+++ b/wk/plot_params_vs_rho.m
@@ -3,6 +3,7 @@ function [ ] = plot_params_vs_rho(geom,prof_folder,rho,rho_min,rho_max,Lref,mref
if FROMPROFILE
% profiles = read_DIIID_profile([prof_folder,'/profiles.txt']);
[params,profiles] = get_param_from_profiles(prof_folder,rho,Lref,mref,Bref,FROMPROFILE);
+ rho_a = profiles.ne.x;
m_e = 9.11e-31;
sigma = sqrt(m_e/mref);
TAU_a = profiles.ti.y./profiles.te.y;
@@ -42,12 +43,15 @@ if FROMPROFILE
else
rho_a = linspace(rho_min,rho_max,100);
- NU_a = zeros(size(rho_a));
- BETA_a = zeros(size(rho_a));
- K_Ne_a = zeros(size(rho_a));
- K_Ti_a = zeros(size(rho_a));
- K_Te_a = zeros(size(rho_a));
-
+ NU_a = zeros(size(rho_a));
+ BETA_a = zeros(size(rho_a));
+ K_Ne_a = zeros(size(rho_a));
+ K_Ti_a = zeros(size(rho_a));
+ K_Te_a = zeros(size(rho_a));
+ Kappa_a = zeros(size(rho_a));
+ Delta_a = zeros(size(rho_a));
+ Zeta_a = zeros(size(rho_a));
+ geom_a = cell(size(rho_a));
for i = 1:numel(rho_a)
rho_ = rho_a(i);
[params,profiles] = get_param_from_profiles(prof_folder,rho_,Lref,mref,Bref,FROMPROFILE);
@@ -56,6 +60,12 @@ else
K_Ne_a(i) = params.K_Ne;
K_Ti_a(i) = params.K_Ti;
K_Te_a(i) = params.K_Te;
+
+ geom_ = get_miller_GENE_py(prof_folder,rho_);
+ Kappa_a(i) = geom_.kappa;
+ Delta_a(i) = geom_.delta;
+ Zeta_a(i) = geom_.zeta;
+ geom_a{i} = geom_;
end
figure
subplot(231)
@@ -88,8 +98,15 @@ else
grid on
end
subplot(133)
- [R,Z] = plot_miller(geom,rho,32,0);
- plot(R,Z,'-b');
+ [R,Z] = plot_miller(geom,rho,128,0);
+ plot(R,Z,'-b'); hold on
+ % rhos = linspace(0.09,0.99,6);
+ % for i = 1:numel(rhos)
+ % rho_ = rhos(i);
+ % geom_ = get_miller_GENE_py(prof_folder,rho_);
+ % [R,Z] = plot_miller(geom_,rho_,128,0);
+ % plot(R,Z,'-k');
+ % end
xlabel('R [m]');
ylabel('Z [m]');
axis tight
diff --git a/wk/plot_spectrum.m b/wk/plot_spectrum.m
new file mode 100644
index 0000000000000000000000000000000000000000..264f3a72b3acb75aaf3e96e22353b17192cf3b6f
--- /dev/null
+++ b/wk/plot_spectrum.m
@@ -0,0 +1,45 @@
+function [FIGURE] = plot_spectrum(data,options)
+
+
+options.INTERP = 0;
+options.POLARPLOT = 0;
+options.AXISEQUAL = 1;
+options.LOGSCALE = 0;
+options.CLIMAUTO = 1;
+options.PLAN = 'kxky'; options.COMP = data.grids.Nz/2+1;
+options.RESOLUTION = 256;
+options.BWR = 0; % bluewhitered plot or gray
+options.COMP = 'avg';
+[toplot] = process_field(data,options);
+
+% Time averaging
+fkykx = squeeze(mean(abs(toplot.FIELD(:,:,:)),3));
+
+FNAME = toplot.FILENAME;
+
+
+kx = toplot.X(1,:);
+Nkx = numel(kx);
+kx = toplot.X(1,Nkx/2:end);
+fkx = squeeze(fkykx(1,Nkx/2:end));
+ky = toplot.Y(:,data.grids.Nkx/2+1);
+fky = squeeze(fkykx(:,data.grids.Nkx/2+1));
+
+FIGURE.fig = figure;
+subplot(121)
+ scale = max(fkx)*options.NORMALIZE + 1*(1-options.NORMALIZE);
+ plot(kx,fkx./scale)
+ xlabel(toplot.XNAME)
+ ylabel(['$|',toplot.FIELDNAME,'|$'])
+
+
+subplot(122)
+ scale = max(fky)*options.NORMALIZE + 1*(1-options.NORMALIZE);
+ plot(ky,fky./scale)
+ xlabel(toplot.YNAME)
+ ylabel(['$|',toplot.FIELDNAME,'|$'])
+
+
+FIGURE.FIGNAME = [FNAME,'_spectrum'];
+
+end
\ No newline at end of file
diff --git a/wk/plot_toroidal_fluxtube_geom.m b/wk/plot_toroidal_fluxtube_geom.m
index a43c212cca0c60b176f70f2f63e1598d201a6ec4..3c2dc8c5a03f573561a98725c34302e2b90a0a42 100644
--- a/wk/plot_toroidal_fluxtube_geom.m
+++ b/wk/plot_toroidal_fluxtube_geom.m
@@ -1,3 +1,8 @@
+%% This script is made for illustrating different magnetic geometry
+% it is not bugfree but work pretty well for illustration purpose of the
+% magnetic flux surface of a tokamak and for the particle trajectory in it.
+% ! The particle trajectory is not correct (does not work for Z-pinch for
+% ! example) but is good enough for illustation.
gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory
addpath(genpath([gyacomodir,'matlab'])) % ... add
addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
@@ -9,25 +14,29 @@ default_plots_options
Nx = 128;
Ny = 128;
Nz = 64;
-Nturns = 1;
-x = linspace(-60,60,Nx)*0.001;
+Nturns = 1.0;
+lr = 0.05; %larmor radius
+Wc = 750; %~cyclotron freq
+vd = 0.25; %~drift velocty
+Gx = linspace(-60,60,Nx)*0.001;
y = linspace(-60,60,Ny)*0.001;
FIELD = ones(Nx,Ny,Nz);
z = linspace(-Nturns*pi,Nturns*pi,Nz);
-N_field_lines = 120;
+N_field_lines = 10;
N_magn_flux_surf = 1;
openangle = pi/3;
-phi0 = openangle/2; phi1= 2*pi-openangle/2;
+% phi0 = openangle/2; phi1= 2*pi-openangle/2;
+phi0 = 0; phi1= 2*pi-1.3;
PLT_FTUBE = 0;
PLT_BASIS = 0;
PLT_DATA = 0;
R0 = 1; % Torus major radius
Z0 = 0;
xpoint = 0;
-select = 3;
+select = 1;
switch select
case 1
- % CBC
+ % Z-pinch
rho = 0.50; drho = 0.1;% Torus minor radius
eps = 0.3;
q0 = 0.000001; % Flux tube safety factor
@@ -78,36 +87,20 @@ switch select
zeta = 0.1;
s_zeta = 0.0;
case 5
- % Fake x-point DIII-D
- N_magn_flux_surf = 2;
- phi0 = 0; phi1= pi/48;
- N_field_lines = 0;
- rho = 0.40; drho = 0.1;% Torus minor radius
- eps = 0.3;
- q0 = 4.8; % Flux tube safety factor
- shear = 2.5;
- kappa = 1.0;
- s_kappa= 0.0;
- delta = 0.0;
- s_delta= 0.0;
- zeta = 0.0;
- s_zeta = 0.0;
- xpoint = 0.5;
- case 6
% play with DIII-D edge
rho = 0.90; drho = 0.1;% Torus minor radius
eps = 0.3;
- q0 = 1.0; % Flux tube safety factor
- Nturns = max(1,q0);
+ q0 = 1.8; % Flux tube safety factor
+ % Nturns = max(1,q0);
shear = 2.5;
- kappa = 0.5;
+ kappa = 1.5;
s_kappa= 1.0;
- delta = 0.0;
+ delta = 0.3;
s_delta= 0.0;
- zeta =-1.0;
+ zeta = 0.0;
s_zeta = 0.0;
end
-
+% Nturns = 1/q0;
Nptor = 64;
% Toroidal angle for the flux surf
% phi = linspace(0+openangle/2, 2*pi-openangle/2, Nptor);
@@ -115,20 +108,23 @@ phi = linspace(phi0,phi1, Nptor);
theta = linspace(-pi, pi, Nptor) ; % Poloidal angle
[p, t] = meshgrid(phi, theta);
Tgeom = 1;
-
DIMENSIONS = [600 600 1200 600];
rho = rho*eps; drho = drho*eps;
+iota = 1/q0;
% field line coordinates
Xfl = @(z) (R0+rho*cos(z+asin(delta*sin(z)))).*cos(q0*z);
Yfl = @(z) (R0+rho*cos(z+asin(delta*sin(z)))).*sin(q0*z);
Zfl = @(z) Z0+kappa*rho*sin(z+zeta*sin(2*z)+xpoint*cos(2*z));
Rvec= @(z) [Xfl(z); Yfl(z); Zfl(z)];
% xvec shearless
-xX = @(z) (Xfl(z)-R0*cos(q0*z))./sqrt((Xfl(z)-R0*cos(q0*z)).^2+(Yfl(z)-R0*sin(q0*z)).^2+Zfl(z).^2);
-xY = @(z) (Yfl(z)-R0*sin(q0*z))./sqrt((Xfl(z)-R0*cos(q0*z)).^2+(Yfl(z)-R0*sin(q0*z)).^2+Zfl(z).^2);
-xZ = @(z) Zfl(z)./sqrt((Xfl(z)-R0*cos(q0*z)).^2+(Yfl(z)-R0*sin(q0*z)).^2+Zfl(z).^2);
+% xX = @(z) (Xfl(z)-R0*cos(q0*z))./sqrt((Xfl(z)-R0*cos(q0*z)).^2+(Yfl(z)-R0*sin(q0*z)).^2+Zfl(z).^2);
+% xY = @(z) (Yfl(z)-R0*sin(q0*z))./sqrt((Xfl(z)-R0*cos(q0*z)).^2+(Yfl(z)-R0*sin(q0*z)).^2+Zfl(z).^2);
+% xZ = @(z) Zfl(z)./sqrt((Xfl(z)-R0*cos(q0*z)).^2+(Yfl(z)-R0*sin(q0*z)).^2+Zfl(z).^2);
+xX = @(z) cos(z+asin(delta*sin(z))).*cos(q0*z)./sqrt(cos(z+asin(delta*sin(z))).^2+kappa^2*sin(z+zeta*sin(2*z)).^2);
+xY = @(z) cos(z+asin(delta*sin(z))).*sin(q0*z)./sqrt(cos(z+asin(delta*sin(z))).^2+kappa^2*sin(z+zeta*sin(2*z)).^2);
+xZ = @(z) kappa*sin(z+zeta*sin(2*z))./sqrt(cos(z+asin(delta*sin(z))).^2+kappa^2*sin(z+zeta*sin(2*z)).^2);
xvec= @(z) [xX(z); xY(z); xZ(z)];
-% bvec
+% bvec shearless
bX = @(z) Tgeom*(rho*cos(z).*cos(q0*z) - q0*Yfl(z))./sqrt(Tgeom*(rho*cos(z).*cos(q0*z) - q0*Yfl(z)).^2+(rho*cos(z).*sin(q0*z) + q0*Xfl(z)).^2+(rho*cos(z)).^2);
bY = @(z) (rho*cos(z).*sin(q0*z) + q0*Xfl(z))./sqrt(Tgeom*(rho*cos(z).*cos(q0*z) - q0*Yfl(z)).^2+(rho*cos(z).*sin(q0*z) + q0*Xfl(z)).^2+(rho*cos(z)).^2);
bZ = @(z) rho*cos(z)./sqrt(Tgeom*(rho*cos(z).*cos(q0*z) - q0*Yfl(z)).^2+(rho*cos(z).*sin(q0*z) + q0*Xfl(z)).^2+(rho*cos(z)).^2);
@@ -141,13 +137,6 @@ yvec= @(z) [yX(z); yY(z); yZ(z)];
scale = 0.10;
-% Plot plane result
-OPTIONS.POLARPLOT = 0;
-OPTIONS.PLAN = 'xy';
-rhostar = 0.1;
-[X,Y] = meshgrid(x,y);
-max_ = 0;
-
figure; set(gcf, 'Position', DIMENSIONS)
%plot magnetic geometry
@@ -171,7 +160,8 @@ figure; set(gcf, 'Position', DIMENSIONS)
end
%plot field lines
if N_field_lines > 0
- theta = linspace(-Nturns*pi, Nturns*pi, 256) ; % Poloidal angle
+ % theta = linspace(-Nturns*pi, Nturns*pi, 1024) ; % Poloidal angle
+ theta = linspace(0, Nturns*2*pi, 1024) ; % Poloidal angle
vecfl = Rvec(theta);
plot3(vecfl(1,:),vecfl(2,:),vecfl(3,:),'-b'); hold on;
% Multiple field lines
@@ -179,17 +169,22 @@ figure; set(gcf, 'Position', DIMENSIONS)
for ifl = 1:N_field_lines-1
% rotation for multiple field lines
t = q0*pi*ifl;
- x_ = [x_ cos(t)*vecfl(1,:) - sin(t)*vecfl(2,:)];
- y_ = [y_ sin(t)*vecfl(1,:) + cos(t)*vecfl(2,:)];
- z_ = [z_ vecfl(3,:)];
+ % x_ = [x_ cos(t)*vecfl(1,:) - sin(t)*vecfl(2,:)];
+ % y_ = [y_ sin(t)*vecfl(1,:) + cos(t)*vecfl(2,:)];
+ % z_ = [z_ vecfl(3,:)];
+ x_ = cos(t)*vecfl(1,:) - sin(t)*vecfl(2,:);
+ y_ = sin(t)*vecfl(1,:) + cos(t)*vecfl(2,:);
+ z_ = vecfl(3,:);
+ % plot3(x_,y_,z_,'-','color',[1.0 0.6 0.6]*0.8); hold on;
+ plot3(x_,y_,z_,'-','color',[0.0 0.5 1.0]*0.8,'LineWidth',1); hold on;
end
- plot3(x_,y_,z_,'.','color',[1.0 0.6 0.6]*0.8); hold on;
+ %plot3(x_,y_,z_,'-','color',[1.0 0.6 0.6]*0.8); hold on;
end
%plot fluxe tube
if PLT_FTUBE
theta = linspace(-Nturns*pi, Nturns*pi, 64) ; % Poloidal angle
%store the shifts in an order (top left to bottom right)
- s_x = rhostar*[x(1) x(end) x(1) x(end)];
+ s_x = rhostar*[Gx(1) Gx(end) Gx(1) Gx(end)];
s_y = rhostar*[y(1) y(1) y(end) y(end)];
for i_ = 1:4
vx_ = Xfl(theta) + s_x(i_)*xX(theta) + s_y(i_)*yX(theta);
@@ -206,19 +201,49 @@ figure; set(gcf, 'Position', DIMENSIONS)
quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*yX(theta),scale*yY(theta),scale*yZ(theta),0,'g');
quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*bX(theta),scale*bY(theta),scale*bZ(theta),0,'b');
end
- xlabel('X');ylabel('Y');zlabel('Z');
- %Plot time dependent data
- if PLT_DATA
- for iz = 1:Nz
- z_ = z(iz);
- Xp = Xfl(z_) + X*xX(z_) + Y*yX(z_);
- Yp = Yfl(z_) + X*xY(z_) + Y*yY(z_);
- Zp = Zfl(z_) + X*xZ(z_) + Y*yZ(z_);
- s=surface(Xp,Yp,Zp,FIELD(:,:,iz)/max(max(max(abs(FIELD)))));
- s.EdgeColor = 'none';
- colormap(bluewhitered);
+ % Particle trajectory
+ if 1
+ % Particle gyration
+ Gx = @(z) lr*cos(z).*cos(Wc*z);
+ Gy = @(z) lr*sin(z).*cos(Wc*z);
+ Gz = @(z) lr*sin(Wc*z);
+ % This is the "real" helicoidal trajectory
+ theta = linspace(0, Nturns*2*pi, 10000) ; % Poloidal angle
+ dtheta = theta(2)-theta(1);
+ dX = 0*theta; dY = dX; dZ = dX; % init drifts
+ for i = 2:numel(theta)
+ th_old = theta(i-1);
+ th_ = theta(i);
+ % find magnetic displacement vector
+ bx_= Xfl(th_)-Xfl(th_old);
+ by_= Yfl(th_)-Yfl(th_old);
+ bz_= Zfl(th_)-Zfl(th_old);
+ % find binormal
+ yx_= xY(th_).*bz_ - xZ(th_).*by_;
+ yy_= xZ(th_).*bx_ - xX(th_).*bz_;
+ yz_= xX(th_).*by_ - xY(th_).*bx_;
+ % normalize
+ yx_=yx_/sqrt(yx_^2+yy_^2+yz_^2);
+ yy_=yy_/sqrt(yx_^2+yy_^2+yz_^2);
+ yz_=yz_/sqrt(yx_^2+yy_^2+yz_^2);
+ % add drift
+ dX(i) = dX(i-1) - vd*yx_*dtheta;
+ dY(i) = dY(i-1) - vd*yy_*dtheta;
+ dZ(i) = dZ(i-1) - vd*yz_*dtheta;
end
+ % add gyration
+ dX = dX + Gx(theta);
+ dY = dY + Gy(theta);
+ dZ = dZ + Gz(theta);
+ % add fieldline
+ Tx = Xfl(theta)+dX;
+ Ty = Yfl(theta)+dY;
+ Tz = Zfl(theta)+dZ;
+ plot3(Tx,Ty,Tz,'-r','LineWidth',1.2); hold on;
+ plot3(Tx(1),Ty(1),Tz(1),'xk','MarkerSize',10); hold on;
+ plot3(Tx(end),Ty(end),Tz(end),'ok','MarkerSize',8,'MarkerFaceColor','k'); hold on;
end
+ xlabel('X');ylabel('Y');zlabel('Z');
%
axis equal
view([1,-2,1])