diff --git a/matlab/load/quick_gene_load.m b/matlab/load/quick_gene_load.m
index fb1077c048d79353f78e53b8ac7094081a5023f7..69f28ca9310e2960bee2a81fad5b47223e7288ed 100644
--- a/matlab/load/quick_gene_load.m
+++ b/matlab/load/quick_gene_load.m
@@ -76,11 +76,16 @@ path = '/home/ahoffman/gene/linear_CBC_results/';
% fname = 'CBC_salpha_nz_24_nv_scan_nw_16_kine.txt';
% fname = 'CBC_miller_nz_24_nv_scan_nw_16_kine.txt';
%---------- CBC
-% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
-% fname = 'CBC_salpha_nx_8_nz_18_nv_18_nw_8_adiabe.txt';
% fname = 'CBC_salpha_nx_8_nz_18_nv_12_nw_8_adiabe.txt';
-% fname = 'CBC_salpha_nx_8_nz_24_nv_60_nw_30_adiabe.txt';
-fname = 'kT_5.3_salpha_nx_8_nz_24_nv_60_nw_30_adiabe.txt';
+% fname = 'CBC_salpha_nx_8_nz_18_nv_18_nw_8_adiabe.txt';
+
+% fname = 'CBC_salpha_nx_8_nz_24_nv_8_nw_4_adiabe.txt';
+% fname = 'CBC_salpha_nx_8_nz_24_nv_16_nw_8_adiabe.txt';
+% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
+fname = 'CBC_salpha_nx_8_nz_24_nv_48_nw_16_adiabe.txt';
+
+% fname = 'kT_5.3_salpha_nx_8_nz_24_nv_60_nw_30_adiabe.txt';
+
% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_kine.txt';
% fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_adiabe.txt';
% fname = 'CBC_miller_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
diff --git a/matlab/plot/create_film.m b/matlab/plot/create_film.m
index e727e49cc8bd3e19f463248cb62747aae233f45b..031ff2c5cb52bfc64fc10922f3461d6768c7ba64 100644
--- a/matlab/plot/create_film.m
+++ b/matlab/plot/create_film.m
@@ -36,7 +36,7 @@ if hmax == hmin
disp('Warning : h = hmin = hmax = const')
else
% Setup figure frame
-fig = figure('Color','white','Position', toplot.DIMENSIONS.*[0 0 1.2 1]);
+fig = figure('Color','white');%,'Position', toplot.DIMENSIONS.*[0.5 0.5 1.0 1]);
if ~strcmp(OPTIONS.PLAN,'sx')
pcolor(X,Y,FIELD(:,:,1)); % to set up
if BWR
@@ -70,7 +70,7 @@ fig = figure('Color','white','Position', toplot.DIMENSIONS.*[0 0 1.2 1]);
if toplot.INTERP
shading interp;
end
- set(pclr, 'edgecolor','none'); pbaspect(toplot.ASPECT);
+ set(pclr, 'edgecolor','none'); %pbaspect(toplot.ASPECT);
if BWR
colormap(bluewhitered)
end
diff --git a/matlab/plot/imagesc_custom.m b/matlab/plot/imagesc_custom.m
index 3b57744e91760e2c17592b2ea526d358ceb38b9a..8183267986feee4dc1851770cf80954a72124e71 100644
--- a/matlab/plot/imagesc_custom.m
+++ b/matlab/plot/imagesc_custom.m
@@ -5,8 +5,8 @@ function [ fig ] = imagesc_custom( XX,YY,FF)
fig = imagesc( XX(1,:), YY(:,1), FF);
set(gca,'YDir','normal')
xticks(XX(1,:));
-% xticklabels(num2str(XX(1,:)));
+xticklabels(num2str(XX(1,:)));
yticks(YY(:,1));
-% yticklabels(num2str(YY(:,1)));
+yticklabels(num2str(YY(:,1)));
end
diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m
index b77806a7bc184ca99056e25e859926aaff614924..a33afa4e2e294e8f2fea3b762254b4a86dd69fd2 100644
--- a/matlab/plot/show_moments_spectrum.m
+++ b/matlab/plot/show_moments_spectrum.m
@@ -70,32 +70,54 @@ for ia = 1:DATA.inputs.Na
end
% plots
% scaling
- if OPTIONS.NORMALIZED
- plt = @(x,i) reshape(x(i,:)./max(x(i,:)),numel(p2j),numel(Time_));
+ if OPTIONS.TAVG_2D
+ nt_half = ceil(numel(Time_)/2);
+ Napjz_tavg = mean(Napjz(:,:,nt_half:end),3);
+ x_ = DATA.grids.Parray;
+ y_ = DATA.grids.Jarray;
+ if OPTIONS.TAVG_2D_CTR
+ [YY,XX] = meshgrid(y_,x_);
+ surf(XX,YY,Napjz_tavg);
+ else
+ imagesc(x_,y_,Napjz_tavg');
+ set(gca,'YDir','normal')
+ xlabel('$p$'); ylabel('$j$');
+ end
+ clb = colorbar; colormap(jet);
+ clb.Label.String = '$\log\langle E(p,j) \rangle_t$';
+ clb.TickLabelInterpreter = 'latex';
+ clb.Label.Interpreter = 'latex';
+ clb.Label.FontSize= 18;
else
- plt = @(x,i) reshape(x(i,:),numel(p2j),numel(Time_));
- end
- if OPTIONS.ST
- imagesc(Time_,p2j,plt(Na_ST,1:numel(p2j)));
- set(gca,'YDir','normal')
- xlabel('$t$'); ylabel(yname);
- if ~isempty(ticks_labels)
- yticks(p2j);
- yticklabels(ticks_labels)
+ if OPTIONS.NORMALIZED
+ plt = @(x,i) reshape(x(i,:)./max(x(i,:)),numel(p2j),numel(Time_));
+ else
+ plt = @(x,i) reshape(x(i,:),numel(p2j),numel(Time_));
end
- if OPTIONS.FILTER
- caxis([0 0.2]);
- title('Relative fluctuation from the average');
+ Nlabelmax = 15;
+ nskip = floor(numel(p2j)/Nlabelmax);
+ if OPTIONS.ST
+ imagesc(Time_,p2j,plt(Na_ST,1:numel(p2j)));
+ set(gca,'YDir','normal')
+ xlabel('$t$'); ylabel(yname);
+ if ~isempty(ticks_labels)
+ yticks(p2j(1:nskip:end));
+ yticklabels(ticks_labels(1:nskip:end))
end
- colorbar
- else
- colors_ = jet(numel(p2j));
- for i = 1:numel(p2j)
- semilogy(Time_,squeeze(Na_ST(i,:)),...
- 'DisplayName',ticks_labels{i},...
- 'color',colors_(i,:)); hold on;
+ if OPTIONS.FILTER
+ caxis([0 0.2]);
+ title('Relative fluctuation from the average');
+ end
+ colorbar
+ else
+ colors_ = jet(numel(p2j));
+ for i = 1:numel(p2j)
+ semilogy(Time_,squeeze(Na_ST(i,:)),...
+ 'DisplayName',ticks_labels{i},...
+ 'color',colors_(i,:)); hold on;
+ end
+ title(plotname)
end
- title(plotname)
end
top_title(DATA.paramshort)
diff --git a/wk/CBC_P_J_scan.m b/wk/CBC_P_J_scan.m
new file mode 100644
index 0000000000000000000000000000000000000000..9d39328fbe612078a356436b8cff412c8a542c0d
--- /dev/null
+++ b/wk/CBC_P_J_scan.m
@@ -0,0 +1,224 @@
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
+EXECNAME = 'gyacomo23_sp';
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%
+SIMID = 'p2_linear_new'; % Name of the simulation
+RERUN = 0; % rerun if the data does not exist
+RUN = 0;
+NU = 1e-3;
+K_T= 6.96;
+P_a = 2:2:40;
+J_a = 2:2:40;
+% collision setting
+CO = 'DG';
+GKCO = 0; % gyrokinetic operator
+COLL_KCUT = 1.75;
+% model
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 1e-4; % electron plasma beta
+% background gradients setting
+K_N = 2.22; % Density '''
+% Geometry
+% GEOMETRY= 'miller';
+GEOMETRY= 's-alpha';
+SHEAR = 0.8; % magnetic shear
+% time and numerical grid
+DT0 = 1e-2;
+TMAX = 30;
+kymin = 0.3;
+NY = 2;
+% arrays for the result
+g_ky = zeros(numel(P_a),numel(J_a),NY/2+1);
+g_avg= g_ky*0;
+g_std= g_ky*0;
+% Naming of the collision operator
+if GKCO
+ CONAME = [CO,'GK'];
+else
+ CONAME = [CO,'DK'];
+end
+
+j = 1;
+for P = P_a
+i = 1;
+for J = J_a
+ %% PHYSICAL PARAMETERS
+ TAU = 1.0; % e/i temperature ratio
+ % SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+ SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+ K_Te = K_T; % ele Temperature '''
+ K_Ti = K_T; % ion Temperature '''
+ K_Ne = K_N; % ele Density '''
+ K_Ni = K_N; % ion Density gradient drive
+ %% GRID PARAMETERS
+ DT = DT0/sqrt(J);
+ PMAX = P; % Hermite basis size
+ JMAX = J; % Laguerre "
+ NX = 8; % real space x-gridpoints
+ LX = 2*pi/0.8; % Size of the squared frequency domain
+ LY = 2*pi/kymin; % Size of the squared frequency domain
+ NZ = 24; % number of perpendicular planes (parallel grid)
+ NPOL = 1;
+ SG = 0; % Staggered z grids option
+ NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+ %% GEOMETRY
+ % GEOMETRY= 's-alpha';
+ EPS = 0.18; % inverse aspect ratio
+ Q0 = 1.4; % safety factor
+ KAPPA = 1.0; % elongation
+ DELTA = 0.0; % triangularity
+ ZETA = 0.0; % squareness
+ PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
+ SHIFT_Y = 0.0;
+ %% TIME PARMETERS
+ DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+ DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+ DTSAVE3D = 2; % Sampling per time unit for 3D arrays
+ DTSAVE5D = 100; % Sampling per time unit for 5D arrays
+ JOB2LOAD = -1; % Start a new simulation serie
+ %% OPTIONS
+ LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+ % Collision operator
+ ABCO = 1; % INTERSPECIES collisions
+ INIT_ZF = 0; ZF_AMP = 0.0;
+ CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
+ NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
+ KERN = 0; % Kernel model (0 : GK)
+ INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
+ NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
+ %% OUTPUTS
+ W_DOUBLE = 0;
+ W_GAMMA = 1; W_HF = 1;
+ W_PHI = 1; W_NA00 = 1;
+ W_DENS = 0; W_TEMP = 1;
+ W_NAPJ = 0; W_SAPJ = 0;
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % unused
+ NA = 1;
+ ADIAB_E = (NA==1);
+ HD_CO = 0.0; % Hyper diffusivity cutoff ratio
+ MU = 0.0; % Hyperdiffusivity coefficient
+ INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
+ MU_X = MU; %
+ MU_Y = MU; %
+ N_HD = 4;
+ HYP_V = 'none';
+ HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+ DMAX = -1;
+ NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+ NMAX = 0;
+ MU_Z = 1.0; %
+ MU_P = 0.0; %
+ MU_J = 0.0; %
+ LAMBDAD = 0.0;
+ NOISE0 = 1.0e-5; % Init noise amplitude
+ BCKGD0 = 0.0; % Init background
+ k_gB = 1.0;
+ k_cB = 1.0;
+ %% RUN
+ setup
+ % naming
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ % check if data exist to run if no data
+ data_ = {};
+ try
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ Ntime = numel(data_.Ts0D);
+ catch
+ data_.outfilenames = [];
+ end
+ if RUN && (RERUN || isempty(data_.outfilenames) || Ntime < 10)
+ % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
+ % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
+ end
+ try
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+ if numel(data_.Ts3D)>10
+ if numel(data_.Ts3D)>10
+ % Load results after trying to run
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+
+ % linear growth rate (adapted for 2D zpinch and fluxtube)
+ options.TRANGE = [0.5 1]*data_.Ts3D(end);
+ options.NPLOTS = 0; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
+ options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
+
+ [~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window
+ [~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ...
+ field = 0;
+ field_t = 0;
+ for ik = 2:NY/2+1
+ field = squeeze(sum(abs(data_.PHI),3)); % take the sum over z
+ field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
+ to_measure = log(field_t(it1:it2));
+ tw = double(data_.Ts3D(it1:it2));
+ % gr = polyfit(tw,to_measure,1);
+ gr = fit(tw,to_measure,'poly1');
+ err= confint(gr);
+ g_ky(i,j,ik) = gr.p1;
+ g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
+ end
+ [gmax, ikmax] = max(g_ky(i,j,:));
+
+ msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
+ end
+ end
+ catch
+ g_ky(i,j,:) = 0;
+ g_std(i,j,:) = 0;
+ end
+ i = i + 1;
+end
+j = j + 1;
+end
+
+if 0
+%% Check time evolution
+figure;
+to_measure = log(field_t);
+plot(data_.Ts3D,to_measure); hold on
+plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
+end
+
+%% take max growth rate among z coordinate
+y_ = g_ky(:,:,2);
+e_ = g_std(:,:,2);
+
+%%
+if(numel(J_a)>1 && numel(P_a)>1)
+%% Save metadata
+ pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
+ jmin = num2str(min(J_a)); jmax = num2str(max(J_a));
+filename = [num2str(NX),'x',num2str(NZ),'_ky_',num2str(kymin),...
+ '_P_',pmin,'_',pmax,...
+ '_J_',jmin,'_',jmax,'_',CONAME,'_',num2str(NU),'_kT_',num2str(K_T),'.mat'];
+metadata.name = filename;
+metadata.kymin = kymin;
+metadata.title = ['$\nu_{',CONAME,'}=$',num2str(NU),'$\kappa_T=$',num2str(K_T),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)];
+metadata.par = [num2str(NX),'x1x',num2str(NZ)];
+metadata.nscan = 2;
+metadata.s2name = '$P$';
+metadata.s2 = P_a;
+metadata.s1name = '$J$';
+metadata.s1 = J_a;
+metadata.dname = '$\gamma c_s/R$';
+metadata.data = y_;
+metadata.err = e_;
+save([SIMDIR,filename],'-struct','metadata');
+disp(['saved in ',SIMDIR,filename]);
+clear metadata tosave
+end
diff --git a/wk/CBC_kT_PJ_scan.m b/wk/CBC_kT_PJ_scan.m
index 0b428d574d725ee41ba9916601bec0616b1c84eb..12086c65eec5ea82c338855913cddad49f4b0572 100644
--- a/wk/CBC_kT_PJ_scan.m
+++ b/wk/CBC_kT_PJ_scan.m
@@ -10,20 +10,20 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
SIMID = 'p2_linear_new'; % Name of the simulation
RERUN = 0; % rerun if the data does not exist
RUN = 1;
-KT_a = [3:0.5:6.5 6.96];
-% KT_a = [5.5 6];
-% P_a = [25];
+% KT_a = [3:0.5:6.5 6.96];
+KT_a = [6.96];
+P_a = [2];
% P_a = [3:1:29];
-P_a = 2:2:10;
+% P_a = 2:2:10;
J_a = floor(P_a/2);
% collision setting
-CO = 'LD';
-NU = 0.1;
-GKCO = 1; % gyrokinetic operator
+CO = 'DG';
+NU = 1e-3;
+GKCO = 0; % gyrokinetic operator
COLL_KCUT = 1.75;
% model
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
-BETA = 1e-4; % electron plasma beta
+BETA = 0; % electron plasma beta
% background gradients setting
K_N = 2.22; % Density '''
% Geometry
@@ -31,7 +31,7 @@ K_N = 2.22; % Density '''
GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
-DT0 = 1e-2;
+DT0 = 5e-3;
TMAX = 30;
kymin = 0.3;
NY = 2;
@@ -134,13 +134,16 @@ for KT = KT_a
data_ = {};
try
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ if numel(data_.Ts3D) < 10
+ data_.outfilenames = [];
+ end
catch
data_.outfilenames = [];
end
if RUN && (RERUN || isempty(data_.outfilenames))
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
- % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
+ % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
end
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
diff --git a/wk/CBC_kT_nu_scan.m b/wk/CBC_kT_nu_scan.m
index d851b9c9e1a5e4c343a46d3841d7118b4e497ebd..c211ed46076f17322b291b548e16ab94d50be021 100644
--- a/wk/CBC_kT_nu_scan.m
+++ b/wk/CBC_kT_nu_scan.m
@@ -9,16 +9,17 @@ EXECNAME = 'gyacomo23_sp';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%
SIMID = 'p2_linear_new'; % Name of the simulation
-RERUN = 1; % rerun if the data does not exist
+RERUN = 0; % rerun if the data does not exist
RUN = 1;
-KT_a = [3.5 4.0 4.5 5.0 5.5 6.0 6.5 6.96];
-NU_a = [0 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2 0.5];
-% KT_a = 3.5;
-% NU_a = 0.5;
-P = 8;
-J = 4;
+KT_a = [2.5 3.0 3.5 4.0 4.5 5.0];
+% KT_a = [3.0 3.5 4.0 4.5];
+NU_a = [1e-3 2e-3 5e-3 1e-2 2e-2 5e-2 1e-1 2e-1 5e-1 1e+0];
+% KT_a = 4.5;
+% NU_a = 0.01;
+P = 16;
+J = P/2;
% collision setting
-CO = 'LD';
+CO = 'DG';
GKCO = 1; % gyrokinetic operator
COLL_KCUT = 1.75;
% model
@@ -31,7 +32,7 @@ K_N = 2.22; % Density '''
GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
-DT = 1e-2;
+DT = 5e-3;
TMAX = 30;
kymin = 0.3;
NY = 2;
@@ -136,8 +137,8 @@ for KT = KT_a
data_.outfilenames = [];
end
if RUN && (RERUN || isempty(data_.outfilenames))
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
- % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
+ % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
end
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
diff --git a/wk/CBC_ky_PJ_scan.m b/wk/CBC_ky_PJ_scan.m
new file mode 100644
index 0000000000000000000000000000000000000000..0d3a80980a69016ea6421f0341d2c6360f150dbe
--- /dev/null
+++ b/wk/CBC_ky_PJ_scan.m
@@ -0,0 +1,222 @@
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
+EXECNAME = 'gyacomo23_sp';
+% EXECNAME = 'gyacomo23_dp';
+% EXECNAME = 'gyacomo23_debug';
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%
+SIMID = 'p2_linear_new'; % Name of the simulation
+RERUN = 0; % rerun if the data does not exist
+RUN = 1;
+K_T = 5.3;
+P_a = [2 4 8 10];
+% P_a = 10;
+ky_a= [0.05:0.05:1.0];
+% ky_a = 0.05;
+% collision setting
+CO = 'LD';
+GKCO = 1; % gyrokinetic operator
+COLL_KCUT = 1.75;
+NU = 1e-2;
+% model
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 1e-4; % electron plasma beta
+% background gradients setting
+K_N = 2.22; % Density '''
+% Geometry
+% GEOMETRY= 'miller';
+GEOMETRY= 's-alpha';
+SHEAR = 0.8; % magnetic shear
+% time and numerical grid
+DT0 = 5e-3;
+TMAX = 50;
+% arrays for the result
+g_ky = zeros(numel(ky_a),numel(P_a),2);
+g_avg= g_ky*0;
+g_std= g_ky*0;
+% Naming of the collision operator
+if GKCO
+ CONAME = [CO,'GK'];
+else
+ CONAME = [CO,'DK'];
+end
+j = 1;
+for P = P_a
+ J = P/2;
+i = 1;
+for ky = ky_a
+ %% PHYSICAL PARAMETERS
+ TAU = 1.0; % e/i temperature ratio
+ % SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+ SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+ K_Te = K_T; % ele Temperature '''
+ K_Ti = K_T; % ion Temperature '''
+ K_Ne = K_N; % ele Density '''
+ K_Ni = K_N; % ion Density gradient drive
+ %% GRID PARAMETERS
+ DT = DT0/sqrt(J);
+ PMAX = P; % Hermite basis size
+ JMAX = P/2; % Laguerre "
+ NX = 8; % real space x-gridpoints
+ NY = 2;
+ LX = 2*pi/0.8; % Size of the squared frequency domain
+ LY = 2*pi/ky; % Size of the squared frequency domain
+ NZ = 24; % number of perpendicular planes (parallel grid)
+ NPOL = 1;
+ SG = 0; % Staggered z grids option
+ NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+ %% GEOMETRY
+ % GEOMETRY= 's-alpha';
+ EPS = 0.18; % inverse aspect ratio
+ Q0 = 1.4; % safety factor
+ KAPPA = 1.0; % elongation
+ DELTA = 0.0; % triangularity
+ ZETA = 0.0; % squareness
+ PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
+ SHIFT_Y = 0.0;
+ %% TIME PARMETERS
+ DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+ DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+ DTSAVE3D = 2; % Sampling per time unit for 3D arrays
+ DTSAVE5D = 100; % Sampling per time unit for 5D arrays
+ JOB2LOAD = -1; % Start a new simulation serie
+ %% OPTIONS
+ LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+ % Collision operator
+ ABCO = 1; % INTERSPECIES collisions
+ INIT_ZF = 0; ZF_AMP = 0.0;
+ CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
+ NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
+ KERN = 0; % Kernel model (0 : GK)
+ INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
+ NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
+ %% OUTPUTS
+ W_DOUBLE = 0;
+ W_GAMMA = 1; W_HF = 1;
+ W_PHI = 1; W_NA00 = 1;
+ W_DENS = 0; W_TEMP = 1;
+ W_NAPJ = 0; W_SAPJ = 0;
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % unused
+ NA = 1;
+ ADIAB_E = (NA==1);
+ HD_CO = 0.0; % Hyper diffusivity cutoff ratio
+ MU = 0.0; % Hyperdiffusivity coefficient
+ INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
+ MU_X = MU; %
+ MU_Y = MU; %
+ N_HD = 4;
+ HYP_V = 'none';
+ HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+ DMAX = -1;
+ NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+ NMAX = 0;
+ MU_Z = 1.0; %
+ MU_P = 0.0; %
+ MU_J = 0.0; %
+ LAMBDAD = 0.0;
+ NOISE0 = 1.0e-5; % Init noise amplitude
+ BCKGD0 = 0.0; % Init background
+ k_gB = 1.0;
+ k_cB = 1.0;
+ %% RUN
+ setup
+ % naming
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ % check if data exist to run if no data
+ data_ = {};
+ try
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ Ntime = numel(data_.Ts0D);
+ catch
+ data_.outfilenames = [];
+ end
+ if RUN && (RERUN || isempty(data_.outfilenames) || Ntime < 10)
+ % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 2 2 1 0; cd ../../../wk'])
+ % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
+ end
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+ if numel(data_.Ts3D)>10
+ if numel(data_.Ts3D)>10
+ % Load results after trying to run
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+
+ % linear growth rate (adapted for 2D zpinch and fluxtube)
+ options.TRANGE = [0.5 1]*data_.Ts3D(end);
+ options.NPLOTS = 0; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
+ options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
+
+ [~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window
+ [~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ...
+ field = 0;
+ field_t = 0;
+ for ik = 2:NY/2+1
+ field = squeeze(sum(abs(data_.PHI),3)); % take the sum over z
+ field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
+ to_measure = log(field_t(it1:it2));
+ tw = double(data_.Ts3D(it1:it2));
+ % gr = polyfit(tw,to_measure,1);
+ gr = fit(tw,to_measure,'poly1');
+ err= confint(gr);
+ g_ky(i,j,ik) = gr.p1;
+ g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
+ end
+ [gmax, ikmax] = max(g_ky(i,j,:));
+
+ msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
+ end
+ end
+
+ i = i + 1;
+end
+j = j + 1;
+end
+
+if 0
+%% Check time evolution
+figure;
+to_measure = log(field_t);
+plot(data_.Ts3D,to_measure); hold on
+plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
+end
+
+%% take max growth rate among z coordinate
+y_ = g_ky(:,:,2);
+e_ = g_std(:,:,2);
+
+%%
+if(numel(ky_a)>1 && numel(P_a)>1)
+%% Save metadata
+ pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
+ kymin = num2str(min(ky_a)); kymax= num2str(max(ky_a));
+filename = [num2str(NX),'x',num2str(NZ),'_ky_',kymin,'_',kymax,...
+ '_P_',pmin,'_',pmax,'_',CONAME,'_',num2str(NU),'_kT_',num2str(K_T),'.mat'];
+metadata.name = filename;
+metadata.kymin = ky;
+metadata.title = ['$\nu_{',CONAME,'}=$',num2str(NU),'$\kappa_T=$',num2str(K_T),', $\kappa_N=$',num2str(K_N)];
+metadata.par = [num2str(NX),'x1x',num2str(NZ)];
+metadata.nscan = 2;
+metadata.s2name = '$P$';
+metadata.s2 = P_a;
+metadata.s1name = '$ky$';
+metadata.s1 = ky_a;
+metadata.dname = '$\gamma c_s/R$';
+metadata.data = y_;
+metadata.err = e_;
+save([SIMDIR,filename],'-struct','metadata');
+disp(['saved in ',SIMDIR,filename]);
+clear metadata tosave
+end
diff --git a/wk/CBC_ky_PJ_scan_colless.m b/wk/CBC_ky_PJ_scan_colless.m
new file mode 100644
index 0000000000000000000000000000000000000000..37779f3fc8797a637df89fc069cd2f888ecdc5af
--- /dev/null
+++ b/wk/CBC_ky_PJ_scan_colless.m
@@ -0,0 +1,221 @@
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
+% EXECNAME = 'gyacomo23_sp';
+EXECNAME = 'gyacomo23_dp';
+% EXECNAME = 'gyacomo23_debug';
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%
+SIMID = 'p2_linear_new'; % Name of the simulation
+RERUN = 0; % rerun if the data does not exist
+RUN = 1;
+K_T = 6.96;
+% P_a = [2 4 8 16 32 48];
+P_a = 64;
+% ky_a= [0.05:0.05:1.0];
+ky_a= [0.25:0.05:0.65];
+% collision setting
+CO = 'DG';
+GKCO = 0; % gyrokinetic operator
+COLL_KCUT = 1.75;
+NU = 1e-3;
+% model
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 1e-4; % electron plasma beta
+% background gradients setting
+K_N = 2.22; % Density '''
+% Geometry
+% GEOMETRY= 'miller';
+GEOMETRY= 's-alpha';
+SHEAR = 0.8; % magnetic shear
+% time and numerical grid
+DT0 = 2e-3;
+TMAX = 30;
+% arrays for the result
+g_ky = zeros(numel(ky_a),numel(P_a),NY/2+1);
+g_avg= g_ky*0;
+g_std= g_ky*0;
+% Naming of the collision operator
+if GKCO
+ CONAME = [CO,'GK'];
+else
+ CONAME = [CO,'DK'];
+end
+j = 1;
+for P = P_a
+i = 1;
+for ky = ky_a
+ %% PHYSICAL PARAMETERS
+ TAU = 1.0; % e/i temperature ratio
+ % SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+ SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+ K_Te = K_T; % ele Temperature '''
+ K_Ti = K_T; % ion Temperature '''
+ K_Ne = K_N; % ele Density '''
+ K_Ni = K_N; % ion Density gradient drive
+ %% GRID PARAMETERS
+ DT = DT0/sqrt(J);
+ PMAX = P; % Hermite basis size
+ JMAX = P/2; % Laguerre "
+ NX = 8; % real space x-gridpoints
+ NY = 2;
+ LX = 2*pi/0.8; % Size of the squared frequency domain
+ LY = 2*pi/ky; % Size of the squared frequency domain
+ NZ = 24; % number of perpendicular planes (parallel grid)
+ NPOL = 1;
+ SG = 0; % Staggered z grids option
+ NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+ %% GEOMETRY
+ % GEOMETRY= 's-alpha';
+ EPS = 0.18; % inverse aspect ratio
+ Q0 = 1.4; % safety factor
+ KAPPA = 1.0; % elongation
+ DELTA = 0.0; % triangularity
+ ZETA = 0.0; % squareness
+ PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
+ SHIFT_Y = 0.0;
+ %% TIME PARMETERS
+ DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+ DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+ DTSAVE3D = 2; % Sampling per time unit for 3D arrays
+ DTSAVE5D = 100; % Sampling per time unit for 5D arrays
+ JOB2LOAD = -1; % Start a new simulation serie
+ %% OPTIONS
+ LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+ % Collision operator
+ ABCO = 1; % INTERSPECIES collisions
+ INIT_ZF = 0; ZF_AMP = 0.0;
+ CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
+ NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
+ KERN = 0; % Kernel model (0 : GK)
+ INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
+ NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
+ %% OUTPUTS
+ W_DOUBLE = 0;
+ W_GAMMA = 1; W_HF = 1;
+ W_PHI = 1; W_NA00 = 1;
+ W_DENS = 0; W_TEMP = 1;
+ W_NAPJ = 0; W_SAPJ = 0;
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % unused
+ NA = 1;
+ ADIAB_E = (NA==1);
+ HD_CO = 0.0; % Hyper diffusivity cutoff ratio
+ MU = 0.0; % Hyperdiffusivity coefficient
+ INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
+ MU_X = MU; %
+ MU_Y = MU; %
+ N_HD = 4;
+ HYP_V = 'none';
+ HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+ DMAX = -1;
+ NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+ NMAX = 0;
+ MU_Z = 1.0; %
+ MU_P = 0.0; %
+ MU_J = 0.0; %
+ LAMBDAD = 0.0;
+ NOISE0 = 1.0e-5; % Init noise amplitude
+ BCKGD0 = 0.0; % Init background
+ k_gB = 1.0;
+ k_cB = 1.0;
+ %% RUN
+ setup
+ % naming
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ % check if data exist to run if no data
+ data_ = {};
+ try
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ Ntime = numel(data_.Ts0D);
+ catch
+ data_.outfilenames = [];
+ end
+ if RUN && (RERUN || isempty(data_.outfilenames) || Ntime < 10)
+ % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
+ % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 2 2 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
+ end
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+ if numel(data_.Ts3D)>10
+ if numel(data_.Ts3D)>10
+ % Load results after trying to run
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+
+ % linear growth rate (adapted for 2D zpinch and fluxtube)
+ options.TRANGE = [0.5 1]*data_.Ts3D(end);
+ options.NPLOTS = 0; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
+ options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
+
+ [~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window
+ [~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ...
+ field = 0;
+ field_t = 0;
+ for ik = 2:NY/2+1
+ field = squeeze(sum(abs(data_.PHI),3)); % take the sum over z
+ field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
+ to_measure = log(field_t(it1:it2));
+ tw = double(data_.Ts3D(it1:it2));
+ % gr = polyfit(tw,to_measure,1);
+ gr = fit(tw,to_measure,'poly1');
+ err= confint(gr);
+ g_ky(i,j,ik) = gr.p1;
+ g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
+ end
+ [gmax, ikmax] = max(g_ky(i,j,:));
+
+ msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
+ end
+ end
+
+ i = i + 1;
+end
+j = j + 1;
+end
+
+if 0
+%% Check time evolution
+figure;
+to_measure = log(field_t);
+plot(data_.Ts3D,to_measure); hold on
+plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
+end
+
+%% take max growth rate among z coordinate
+y_ = g_ky(:,:,2);
+e_ = g_std(:,:,2);
+
+%%
+if(numel(ky_a)>1 && numel(P_a)>1)
+%% Save metadata
+ pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
+ kymin = num2str(min(ky_a)); kymax= num2str(max(ky_a));
+filename = [num2str(NX),'x',num2str(NZ),'_ky_',kymin,'_',kymax,...
+ '_P_',pmin,'_',pmax,'_',CONAME,'_',num2str(NU),'_kT_',num2str(K_T),'.mat'];
+metadata.name = filename;
+metadata.kymin = ky;
+metadata.title = ['$\nu_{',CONAME,'}=$',num2str(NU),'$\kappa_T=$',num2str(K_T),', $\kappa_N=$',num2str(K_N)];
+metadata.par = [num2str(NX),'x1x',num2str(NZ)];
+metadata.nscan = 2;
+metadata.s2name = '$P$';
+metadata.s2 = P_a;
+metadata.s1name = '$ky$';
+metadata.s1 = ky_a;
+metadata.dname = '$\gamma c_s/R$';
+metadata.data = y_;
+metadata.err = e_;
+save([SIMDIR,filename],'-struct','metadata');
+disp(['saved in ',SIMDIR,filename]);
+clear metadata tosave
+end
diff --git a/wk/CBC_nu_PJ_scan.m b/wk/CBC_nu_PJ_scan.m
new file mode 100644
index 0000000000000000000000000000000000000000..608ca5d36a4c9804e74f9775243eec4a7c104138
--- /dev/null
+++ b/wk/CBC_nu_PJ_scan.m
@@ -0,0 +1,223 @@
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
+EXECNAME = 'gyacomo23_sp';
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%
+SIMID = 'p2_linear_new'; % Name of the simulation
+RERUN = 0; % rerun if the data does not exist
+RUN = 1;
+nu_a = [1e-3,2e-3,5e-3,1e-2,2e-2,5e-2,1e-1];
+% KT_a = [5.5 6];
+% P_a = [25];
+% P_a = [3:1:29];
+P_a = 2:2:4;
+J_a = floor(P_a/2);
+% collision setting
+CO = 'DG';
+K_T = 6.96;
+GKCO = 1; % gyrokinetic operator
+COLL_KCUT = 1.75;
+% model
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 1e-4; % electron plasma beta
+% background gradients setting
+K_N = 2.22; % Density '''
+% Geometry
+% GEOMETRY= 'miller';
+GEOMETRY= 's-alpha';
+SHEAR = 0.8; % magnetic shear
+% time and numerical grid
+DT0 = 1e-2;
+TMAX = 30;
+kymin = 0.3;
+NY = 2;
+% arrays for the result
+g_ky = zeros(numel(nu_a),numel(P_a),NY/2+1);
+g_avg= g_ky*0;
+g_std= g_ky*0;
+% Naming of the collision operator
+if GKCO
+ CONAME = [CO,'GK'];
+else
+ CONAME = [CO,'DK'];
+end
+
+j = 1;
+for P = P_a
+i = 1;
+for NU = nu_a
+ %% PHYSICAL PARAMETERS
+ TAU = 1.0; % e/i temperature ratio
+ % SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+ SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+ K_Te = K_T; % ele Temperature '''
+ K_Ti = K_T; % ion Temperature '''
+ K_Ne = K_N; % ele Density '''
+ K_Ni = K_N; % ion Density gradient drive
+ %% GRID PARAMETERS
+ J = floor(P/2);
+ DT = DT0/sqrt(J);
+ PMAX = P; % Hermite basis size
+ JMAX = J; % Laguerre "
+ NX = 8; % real space x-gridpoints
+ LX = 2*pi/0.8; % Size of the squared frequency domain
+ LY = 2*pi/kymin; % Size of the squared frequency domain
+ NZ = 24; % number of perpendicular planes (parallel grid)
+ NPOL = 1;
+ SG = 0; % Staggered z grids option
+ NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+ %% GEOMETRY
+ % GEOMETRY= 's-alpha';
+ EPS = 0.18; % inverse aspect ratio
+ Q0 = 1.4; % safety factor
+ KAPPA = 1.0; % elongation
+ DELTA = 0.0; % triangularity
+ ZETA = 0.0; % squareness
+ PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
+ SHIFT_Y = 0.0;
+ %% TIME PARMETERS
+ DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+ DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+ DTSAVE3D = 2; % Sampling per time unit for 3D arrays
+ DTSAVE5D = 100; % Sampling per time unit for 5D arrays
+ JOB2LOAD = -1; % Start a new simulation serie
+ %% OPTIONS
+ LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+ % Collision operator
+ ABCO = 1; % INTERSPECIES collisions
+ INIT_ZF = 0; ZF_AMP = 0.0;
+ CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
+ NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
+ KERN = 0; % Kernel model (0 : GK)
+ INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
+ NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
+ %% OUTPUTS
+ W_DOUBLE = 0;
+ W_GAMMA = 1; W_HF = 1;
+ W_PHI = 1; W_NA00 = 1;
+ W_DENS = 0; W_TEMP = 1;
+ W_NAPJ = 0; W_SAPJ = 0;
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % unused
+ ADIAB_E = (NA==1);
+ HD_CO = 0.0; % Hyper diffusivity cutoff ratio
+ MU = 0.0; % Hyperdiffusivity coefficient
+ INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
+ MU_X = MU; %
+ MU_Y = MU; %
+ N_HD = 4;
+ HYP_V = 'none';
+ HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+ DMAX = -1;
+ NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+ NMAX = 0;
+ MU_Z = 1.0; %
+ MU_P = 0.0; %
+ MU_J = 0.0; %
+ LAMBDAD = 0.0;
+ NOISE0 = 1.0e-5; % Init noise amplitude
+ BCKGD0 = 0.0; % Init background
+ k_gB = 1.0;
+ k_cB = 1.0;
+ %% RUN
+ setup
+ % naming
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ % check if data exist to run if no data
+ data_ = {};
+ try
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ catch
+ data_.outfilenames = [];
+ end
+ if RUN && (RERUN || isempty(data_.outfilenames))
+ % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
+ % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
+ end
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+ if numel(data_.Ts3D)>10
+ if numel(data_.Ts3D)>10
+ % Load results after trying to run
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+
+ % linear growth rate (adapted for 2D zpinch and fluxtube)
+ options.TRANGE = [0.5 1]*data_.Ts3D(end);
+ options.NPLOTS = 0; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
+ options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
+
+ [~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window
+ [~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ...
+ field = 0;
+ field_t = 0;
+ for ik = 2:NY/2+1
+ field = squeeze(sum(abs(data_.PHI),3)); % take the sum over z
+ field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
+ to_measure = log(field_t(it1:it2));
+ tw = double(data_.Ts3D(it1:it2));
+ % gr = polyfit(tw,to_measure,1);
+ gr = fit(tw,to_measure,'poly1');
+ err= confint(gr);
+ g_ky(i,j,ik) = gr.p1;
+ g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
+ end
+ [gmax, ikmax] = max(g_ky(i,j,:));
+
+ msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
+ end
+ end
+
+ i = i + 1;
+end
+j = j + 1;
+end
+
+if 0
+%% Check time evolution
+figure;
+to_measure = log(field_t);
+plot(data_.Ts3D,to_measure); hold on
+plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
+end
+
+%% take max growth rate among z coordinate
+y_ = g_ky(:,:,2);
+e_ = g_std(:,:,2);
+
+%%
+if(numel(nu_a)>1 && numel(P_a)>1)
+%% Save metadata
+numin = num2str(min(nu_a)); numax = num2str(max(nu_a));
+ pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
+filename = [num2str(NX),'x',num2str(NZ),'_ky_',num2str(kymin),...
+ '_nu_',numin,'_',numax,...
+ '_P_',pmin,'_',pmax,'_',CONAME,'_kT_',num2str(K_T),'.mat'];
+metadata.name = filename;
+metadata.kymin = kymin;
+metadata.title = ['$\kappa_T=$',num2str(K_T),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)];
+metadata.par = [num2str(NX),'x1x',num2str(NZ)];
+metadata.nscan = 2;
+metadata.s1name = ['$\nu_{',CONAME,'}=$'];
+metadata.s1 = nu_a;
+metadata.s2name = '$P$, $J=\lfloor P/2 \rfloor$';
+metadata.s2 = P_a;
+metadata.dname = '$\gamma c_s/R$';
+metadata.data = y_;
+metadata.err = e_;
+metadata.date = date;
+save([SIMDIR,filename],'-struct','metadata');
+disp(['saved in ',SIMDIR,filename]);
+clear metadata tosave
+end
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 66ce890e7b80045655367f18762bdee987a7ea6b..c36886926c94ac1ba6ad79bcc3df526042a34582 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -53,13 +53,14 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x32x16_Nexc_5_00/';
% folder = '/misc/gene_results/CBC/KT_4.5_192x96x24x30x16_00/';
-folder = '/misc/gene_results/CBC/KT_5.0_192x96x24x30x16_00/';
+% folder = '/misc/gene_results/CBC/KT_5.0_192x96x24x30x16_00/';
% folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x16x8_Nexc_5/';
% folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x8x4_Nexc_5/';
% folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x8x4_Nexc_5_smallvbox/';
% folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x16x8_Nexc_5_largexbox/';
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x16x8_Muz_0.02/';
-
+% folder = '/misc/gene_results/kT_scan_nu0/30x16x128x64x24/kT_4.0/';
+folder = '/misc/gene_results/kT_scan_nu0/42x24x128x64x24/kT_6.5/';
% Miller CBC
% folder = '/misc/gene_results/CBC/Miller_KT_6.96_128x64x24x16x8/';
% folder = '/misc/gene_results/CBC/Miller_KT_4.1_128x64x24x16x8/';
diff --git a/wk/benchmark_and_scan_scripts/Dimits_2000_fig3_dashed_low.txt b/wk/benchmark_and_scan_scripts/Dimits_2000_fig3_dashed_low.txt
new file mode 100644
index 0000000000000000000000000000000000000000..b12b9e67cea955997f9581ad4ede36287d9ca15f
--- /dev/null
+++ b/wk/benchmark_and_scan_scripts/Dimits_2000_fig3_dashed_low.txt
@@ -0,0 +1,24 @@
+6.074498567335244, 0.0664424234873966
+6.275071633237822, 0.6469830987085263
+6.50429799426934, 1.1655306198722162
+6.762177650429798, 1.6634928544732954
+7.048710601719197, 2.2443894712488195
+7.335243553008594, 2.7631742867821103
+7.679083094555873, 3.4064199991694704
+8.13753581661891, 3.9880284990537884
+8.538681948424069, 4.569399704568509
+9.083094555873924, 5.254883814744286
+9.570200573065902, 5.712387359328931
+9.999999999999998, 6.148949675796567
+10.716332378223493, 6.752329934091489
+11.20343839541547, 7.16842561118131
+11.518624641833808, 7.376770067688219
+12.177650429799423, 7.814281561634245
+13.15186246418338, 8.377321777097535
+13.810888252148997, 8.752721469801324
+14.84240687679083, 9.191775377149739
+15.702005730659023, 9.54730166639971
+16.446991404011456, 9.798833698173425
+17.335243553008596, 10.133774700860785
+18.28080229226361, 10.406841196675504
+19.16905444126074, 10.617558596878393
diff --git a/wk/benchmark_and_scan_scripts/Dimits_2000_fig3_full_no_GF.txt b/wk/benchmark_and_scan_scripts/Dimits_2000_fig3_full_no_GF.txt
index 055da52cb7ea0464cb063fada7a1f28672d9160b..7877a59e27f0cf3934d5deb0bfad9c48cc350ef4 100644
--- a/wk/benchmark_and_scan_scripts/Dimits_2000_fig3_full_no_GF.txt
+++ b/wk/benchmark_and_scan_scripts/Dimits_2000_fig3_full_no_GF.txt
@@ -1,21 +1,11 @@
-5.297267537770027 -0.051044540400214444
-4.754512410510305 0.004436559014473929
-5.089866441088001 0.2660247587319642
-5.451261144141954 0.43432466047664775
-5.709658557756645 0.4345657778143899
-6.046499478583758 0.005642145703188106
-5.2693381128147845 0.9193804088144457
-5.914809225953066 1.1625954673959136
-6.948639997749572 1.0515850450989888
-6.327280618385599 1.6108808217278856
-6.791552052209941 2.003226953703461
-6.894710086541032 2.0966358103451412
-6.9714657723892515 2.4512953024314665
-9.039247874651133 2.1732870120136702
-9.036515211490045 3.4423357840231628
-8.905066076197096 4.487314214067991
-8.955620344677234 5.00991193189239
-12.023959026127072 8.073428266587367
-13.934050389504964 9.026999113893783
-16.026145156655943 9.458189248979773
-17.988719727215916 10.038558680927418
+5.3, 0.03534795644626776
+5.5, 0.4494240176739783
+6.0, 0.0578349376676659
+6.0, 1.2154805112829408
+6.5, 1.568171058860658
+6.96, 2.0654884014517894
+6.96, 0.9906106990689576
+6.96, 2.4583399084740396
+9.0, 4.489742780495503
+9.0, 5.027378885908156
+9.0, 2.216032823102413
diff --git a/wk/benchmark_and_scan_scripts/Lin_1999_fig2.txt b/wk/benchmark_and_scan_scripts/Lin_1999_fig2.txt
index 964d1f9ad85c2868eedfd0297c72607f9ab24aac..2f3885b6e7cf5122bb9d7d54eb063c5845829fed 100644
--- a/wk/benchmark_and_scan_scripts/Lin_1999_fig2.txt
+++ b/wk/benchmark_and_scan_scripts/Lin_1999_fig2.txt
@@ -1,4 +1,4 @@
-0, 0.022113821138211365
+0.001, 0.022113821138211365
0.04541040549724383, 0.17300813008130078
0.08992977422034283, 0.2809756097560976
0.1803843539983388, 0.45788617886178873
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index 3da0a52b25139174af76c1474110c671296f3392..e73ed755311964f45701a830945a5f9ece53b3f6 100644
--- a/wk/fast_analysis.m
+++ b/wk/fast_analysis.m
@@ -48,10 +48,19 @@ PARTITION = '/misc/gyacomo23_outputs/';
% resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_CBC/9x2x64x48x16/nu_0.1';
% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/9x5x128x64x24/nu_0.5';
% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/5x3x128x64x24/nu_0.5';
-resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x128x64x24/nu_0.5';
+% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x128x64x24/nu_0.5';
% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x128x64x24_Lx200/nu_0.5';
-% resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_kT_5.3/9x2x128x64x24/nu_0.01';
+% resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_kT_5.3/9x5x128x64x24/nu_0.005';
+% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/11x6x192x96x24/nu_0.1';
+% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/9x5x192x96x24/nu_0.2';
+% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/17x9x192x96x24/nu_0.2';
+
+% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/9x5x192x96x24/nu_0.02';
+resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/17x9x192x96x24/nu_0.02';
+
+% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x192x96x24/nu_0.2';
+% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/17x9x192x96x24/nu_0.02';
%% kT eff study
% resdir = 'paper_2_GYAC23/kT_eff_study/1x3x128x64x24_kT_3.0/Lx120';
% resdir = 'paper_2_GYAC23/kT_eff_study/1x3x128x64x24_kT_3.0/Lx240';
@@ -61,6 +70,16 @@ resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x128x64x24/nu_0.5
% resdir = 'paper_2_GYAC23/kT_eff_study/7x3x128x64x24_kT_3.6/truncation';
% resdir = 'paper_2_GYAC23/kT_eff_study/7x3x128x64x24_kT_3.6/dmax';
+%% dev
+% PARTITION='';
+% resdir = '/home/ahoffman/gyacomo/testcases/zpinch_3D';
+% resdir = '/home/ahoffman/gyacomo';
+%% CBC benchmark
+% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/3x2x128x64x24/kT_7.0';
+% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/5x3x128x64x24/kT_7.0';
+% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/9x5x128x64x24/kT_7.0';
+% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/17x9x128x64x24/kT_7.0';
+% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/31x16x128x64x24/kT_7.0';
%%
J0 = 00; J1 = 10;
@@ -87,13 +106,15 @@ end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
-[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
+% [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.NAME = '\phi';
+% options.NAME = '\phi';
% options.NAME = '\phi^{NZ}';
% options.NAME = '\omega_z';
-% options.NAME = 'N_i^{00}';
+options.NAME = 'N_i^{00}';
% options.NAME = 's_{Ey}';
% options.NAME = 'n_i^{NZ}';
% options.NAME = 'Q_x';
@@ -104,7 +125,7 @@ options.PLAN = 'xy';
% options.PLAN = 'sx';
options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
-options.TIME = data.Ts3D;
+options.TIME = data.Ts3D(1:1:end);
% options.TIME = [0:1500];
data.EPS = 0.1;
data.a = data.EPS * 2000;
@@ -116,6 +137,7 @@ if 0
%% fields snapshots
% Options
[data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00');
+[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
options.INTERP = 0;
@@ -126,7 +148,7 @@ options.NAME = 'N_i^{00}';
% options.NAME = '\phi';
options.PLAN = 'xy';
options.COMP = 'avg';
-options.TIME = [10 50 80];
+options.TIME = [100 200 300];
options.RESOLUTION = 256;
fig = photomaton(data,options);
% save_figure(data,fig)
@@ -144,10 +166,12 @@ options.ST = 1;
options.NORMALIZED = 0;
options.LOGSCALE = 1;
options.FILTER = 0; %filter the 50% time-average of the spectrum from
+options.TAVG_2D = 1; %Show a 2D plot of the modes, 50% time averaged
+options.TAVG_2D_CTR= 1; %make it contour plot
fig = show_moments_spectrum(data,options);
end
-if 1
+if 0
%% Mode evolution
[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
diff --git a/wk/heat_flux_convergence_analysis.m b/wk/heat_flux_convergence_analysis.m
index 805322675305a2ea226cafeccf90739c06ae2b02..50b6c0f4718a1f3288159668c23e7d73a554cf6b 100644
--- a/wk/heat_flux_convergence_analysis.m
+++ b/wk/heat_flux_convergence_analysis.m
@@ -1,3 +1,5 @@
+figure
+hold on
if 0
%% GYAC Local low res
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/nu=0.05/3x2x64x48x16/CBC/';
@@ -22,7 +24,7 @@ if 0
legend('show')
disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
end
-if 1
+if 0
%% Manually gathered data
QvsPJ = [...
%vp mu Qxav Qxer kT
@@ -36,7 +38,6 @@ QvsPJ = [...
11 06 27.47 05.48 6.96;...
13 05 16.45 04.63 6.96;...
];
-figure
errorbar(QvsPJ(:,1).*QvsPJ(:,2),QvsPJ(:,3),QvsPJ(:,4),'--s',...
'DisplayName','GYAC 64x48x16 ($\nu_{DGDK}=0.001$)'); hold on
end
@@ -62,7 +63,9 @@ if 0
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/7x4x128x64x24/';
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/7x4x128x64x24/CBC_L120/';
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/local_runs/7x4x128x64x24/CBC_L180/';
-% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/21x6x192x96x24/';
+ % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/21x6x192x96x24/';
+ % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/21x6x192x96x24/';
+ % resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/CBC/21x6x192x96x24/';
data = {};
data = compile_results_low_mem(data,resdir,00,10);
% fast heat flux analysis
diff --git a/wk/lin_ITG_salpha.m b/wk/lin_ITG_salpha.m
index 0df72f58186e4edf3ab3930a613f0e299e707dd8..c603b379acbdeb9be97143e08806feabba42b9cc 100644
--- a/wk/lin_ITG_salpha.m
+++ b/wk/lin_ITG_salpha.m
@@ -31,8 +31,8 @@ NA = 1; % number of kinetic species
ADIAB_E = (NA==1); % adiabatic electron model
BETA = 0.0; % electron plasma beta
%% Set up grid parameters
-P = 12;
-J = 6;%P/2;
+P = 2;
+J = 1;%P/2;
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size
NX = 8; % real space x-gridpoints
diff --git a/wk/load_metadata_scan.m b/wk/load_metadata_scan.m
index 728438f4bc39f0e792f5bed7d4e02ad322af50fd..040038705b053b23e0832d9e2b4a01c5c9b5ef8b 100644
--- a/wk/load_metadata_scan.m
+++ b/wk/load_metadata_scan.m
@@ -49,10 +49,27 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_10_LDGK_0.05.mat';
% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_20_DGDK_0.1.mat';
-datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_20_DGGK_0.1.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_20_DGGK_0.1.mat';
% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_20_SGGK_0.1.mat';
% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_10_LDGK_0.1.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_2_J_1_kT_3_5_nu_0.001_0.1_DGGK.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_4_J_2_kT_2.5_4.5_nu_0.001_0.1_DGGK.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_6_J_3_kT_2.5_4.5_nu_0.001_0.1_DGGK.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_8_J_4_kT_2.5_4.5_nu_0.001_0.1_DGGK.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.3_P_16_J_8_kT_2.5_4.5_nu_0.001_0.1_DGGK.mat';
+
+% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_64_DGDK_0.001_kT_6.96.mat';
+%% ky pj scan
+% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_64_DGDK_0.001_kT_6.96.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_48_DGDK_0.001_kT_6.96.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_16_DGGK_0.01_kT_5.3.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_16_SGGK_0.01_kT_5.3.mat';
+datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_10_LDGK_0.01_kT_5.3.mat';
+%% Data for the paper :
+% Collisionless kT threshold
+% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_30_DGDK_0.001.mat';
+
%% Chose if we filter gamma>0.05
FILTERGAMMA = 1;
@@ -63,18 +80,20 @@ if FILTERGAMMA
d.data = d.data.*(d.data>0.025);
d.err = d.err.*(d.data>0.025);
end
-if 1
+if 0
%% Pcolor of the peak
figure;
% [XX_,YY_] = meshgrid(d.s1,d.s2);
[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
pclr=imagesc_custom(XX_,YY_,d.data'.*(d.data>0)');
+% pclr=contourf(1:numel(d.s1),1:numel(d.s2),d.data'.*(d.data>0)');
+% pclr=surf(1:numel(d.s1),1:numel(d.s2),d.data'.*(d.data>0)');
title(d.title);
xlabel(d.s1name); ylabel(d.s2name);
-set(gca,'XTicklabel',d.s1)
-set(gca,'YTicklabel',d.s2)
+set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
+set(gca,'YTick',1:numel(d.s2),'YTicklabel',d.s2)
colormap(jet)
-% colormap(bluewhitered)
+colormap(bluewhitered)
clb=colorbar;
clb.Label.String = '$\gamma c_s/R$';
clb.Label.Interpreter = 'latex';
@@ -99,9 +118,11 @@ xlabel(d.s1name); ylabel(d.dname);title(d.title);
xlim([d.s1(1) d.s1(end)]);
colormap(colors_);
clb = colorbar;
-caxis([d.s2(1)-0.5,d.s2(end)+0.5]);
+% caxis([d.s2(1)-0.5,d.s2(end)+0.5]);
+clim([1 numel(d.s2)+1]);
clb.Ticks=linspace(d.s2(1),d.s2(end),numel(d.s2));
-clb.YTick=d.s2;
+clb.Ticks =1.5:numel(d.s2)+1.5;
+clb.TickLabels=d.s2;
clb.Label.String = d.s2name;
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
@@ -169,29 +190,38 @@ if 0
figure; i_ = 0;
% target_ = 2.72724991618068013377e-01; % Value for nuDGDK = 1.0, kT=6.96, (40,20), Nkx=8
% target_ = 2.79666916212537142172e-01; % Value for nuDGDK = 0.05, kT=6.96, (40,20), Nkx=8
-target_ = 2.73048910051283844069e-01; % Value for nuDGDK = 0.0, kT=6.96, (40,20), Nkx=8
+% target_ = 2.71455e-01; % Value for nuDGDK = 0.001, kT=6.96, (40,20), Nkx=8
+% target_ = 1.39427e-01; % Value for nuDGDK = 0.001, kT=5.3, (30,16), Nkx=8
+target_ = 2.72510405826983714839e-01 % Value for nuDGDK = 0.001, kT=6.96, (50,25), Nkx=8
+% target_ = 2.73048910051283844069e-01; % Value for nuDGDK = 0.0, kT=6.96, (40,20), Nkx=8
% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
% eps_ = log(abs(target_ - d.data)/abs(target_));
eps_ = max(-10,log(abs(target_ - d.data)/abs(target_)));
sign_ = 1;%sign(d.data - target_);
eps_ = d.data;
for i = 1:numel(d.s1)
- target_ = d.data(i,end);
- eps_(i,:) = log(abs(target_ - d.data(i,1:end)));
- if target_ > 0
-% eps_(i,:) = max(-12,log(abs(target_ - d.data(i,1:end))/abs(target_)));
-% eps_(i,:) = log(abs(target_ - d.data(i,1:end))/abs(target_));
-% eps_(i,:) = min(100,100*abs(target_ - d.data(i,1:end))/abs(target_));
- else
+ for j = 1:numel(d.s2)
+ % target_ = d.data(i,end);
+ % target_ = d.data(i,end);
+ % target_ = d.data(end,j);
+ eps_(i,j) = log(abs(target_ - d.data(i,j))/target);
+ if target_ > 0
+ % eps_(i,:) = max(-12,log(abs(target_ - d.data(i,1:end))/abs(target_)));
+ % eps_(i,:) = log(abs(target_ - d.data(i,1:end))/abs(target_));
+ % eps_(i,:) = min(100,100*abs(target_ - d.data(i,1:end))/abs(target_));
+ else
+ end
end
end
[XX_,YY_] = meshgrid(d.s1,d.s2);
[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
pclr=imagesc_custom(XX_,YY_,eps_'.*(d.data>0)'.*sign_');
+% pclr=contourf(1:numel(d.s1),1:numel(d.s2),eps_'.*(d.data>0)'.*sign_',5);
+% pclr=surf(1:numel(d.s1),1:numel(d.s2),eps_'.*(d.data>0)'.*sign_');
title(d.title);
xlabel(d.s1name); ylabel(d.s2name);
-set(gca,'XTicklabel',d.s1)
-set(gca,'YTicklabel',d.s2)
+set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
+set(gca,'YTick',1:numel(d.s2),'YTicklabel',d.s2)
% colormap(jet)
colormap(bluewhitered)
% caxis([-10, 0]);
diff --git a/wk/local_run.m b/wk/local_run.m
index 5ba805aa00da227eef1cbd9376004d71d5ae12ed..1cd20ee94ef4cf234fdc7fb1c3104e6e8648c1a5 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -13,14 +13,14 @@ addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
%% Set simulation parameters
SIMID = 'dbg'; % Name of the simulation
-RUN = 1; % To run or just to load
+RUN = 0; % To run or just to load
default_plots_options
EXECNAME = 'gyacomo23_sp'; % single precision
% EXECNAME = 'gyacomo23_dp'; % double precision
%% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
-NU = 0.05; % Collision frequency
+NU = 0.001; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_Ne = 0*2.22; % ele Density
K_Te = 0*6.96; % ele Temperature
@@ -31,13 +31,13 @@ NA = 1; % number of kinetic species
ADIAB_E = (NA==1); % adiabatic electron model
BETA = 0.0; % electron plasma beta
%% Set up grid parameters
-P = 16;
+P = 60;
J = P/2;
-DT = 1e-2; % Time step
+DT = 1e-2; % Time step
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size
NX = 8; % real space x-gridpoints
-NY = 12; % real space y-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.1; % Size of the squared frequency domain in y direction
NZ = 24; % number of perpendicular planes (parallel grid)
@@ -71,7 +71,8 @@ CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:S
GKCO = 1; % Gyrokinetic operator
ABCO = 1; % INTERSPECIES collisions
INIT_ZF = 0; % Initialize zero-field quantities
-HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+% HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+HRCY_CLOS = 'monomial'; % Closure model for higher order moments
DMAX = -1;
NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
NMAX = 0;
diff --git a/wk/new_CBC_convergence_GENE_GYAC_comparison.m b/wk/new_CBC_convergence_GENE_GYAC_comparison.m
new file mode 100644
index 0000000000000000000000000000000000000000..8723bdd6ce3ab0be34087d1f19544d54157f7af4
--- /dev/null
+++ b/wk/new_CBC_convergence_GENE_GYAC_comparison.m
@@ -0,0 +1,181 @@
+figure
+hold on
+TRTW = [0.5 1];
+ERRBAR = 1;
+%%%%%%%%%%%%%%%%%%%%%%%%
+%% GYAC Marconi standard res
+resdirs = {...
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/3x2x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/5x3x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/save/kT_scan_1e-3/7x4x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/9x5x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/save/kT_scan_1e-3/11x6x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/17x9x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/31x16x128x64x24/kT_7.0/'
+ };
+QvsPJ = zeros(numel(resdirs),4);
+N = numel(resdirs);
+for i = 1:N
+ data = {};
+ data = compile_results_low_mem(data,resdirs{i},00,10);
+ % fast heat flux analysis
+ T = data.Ts0D;
+ Qx = data.HFLUX_X;
+ Trange = T(end)*TRTW;
+ [~,it0] = min(abs(Trange(1)-T));
+ [~,it1] = min(abs(Trange(2)-T));
+ Qavg = mean(Qx(it0:it1));
+ Qstd = std(Qx(it0:it1));
+ disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
+ QvsPJ(i,1) = data.grids.Np;
+ QvsPJ(i,2) = data.grids.Nj;
+ QvsPJ(i,3) = Qavg;
+ QvsPJ(i,4) = Qstd;
+end
+if ERRBAR
+ errorbar(QvsPJ(:,1).*QvsPJ(:,2),QvsPJ(:,3),QvsPJ(:,4),'LineStyle','-',...
+ 'DisplayName','GM, $\nu_{DGDK}=0.001$','Marker','o',...
+ 'MarkerSize',7,'LineWidth',2);
+else
+ plot((QvsPJ(:,1).*QvsPJ(:,2)),QvsPJ(:,3),'-o','DisplayName','GM, $\nu_{DGDK}=0.001$')
+end
+
+%% GYAC Marconi nu = 0.01
+resdirs = {...
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-2/5x3x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-2/7x4x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-2/9x5x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-2/11x6x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-2/13x7x128x64x24/kT_7.0/'
+ };
+QvsPJ = zeros(numel(resdirs),4);
+N = numel(resdirs);
+for i = 1:N
+ data = {};
+ data = compile_results_low_mem(data,resdirs{i},00,10);
+ % fast heat flux analysis
+ T = data.Ts0D;
+ Qx = data.HFLUX_X;
+ Trange = T(end)*TRTW;
+ [~,it0] = min(abs(Trange(1)-T));
+ [~,it1] = min(abs(Trange(2)-T));
+ Qavg = mean(Qx(it0:it1));
+ Qstd = std(Qx(it0:it1));
+ disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
+ QvsPJ(i,1) = data.grids.Np;
+ QvsPJ(i,2) = data.grids.Nj;
+ QvsPJ(i,3) = Qavg;
+ QvsPJ(i,4) = Qstd;
+end
+if ERRBAR
+ errorbar(QvsPJ(:,1).*QvsPJ(:,2),QvsPJ(:,3),QvsPJ(:,4),'LineStyle','none',...
+ 'DisplayName','GM, $\nu_{DGDK}=0.01$','Marker','o',...
+ 'MarkerSize',7,'LineWidth',2);
+else
+ plot((QvsPJ(:,1).*QvsPJ(:,2)),QvsPJ(:,3),'-o','DisplayName','GM, $\nu_{DGDK}=0.01$')
+end
+
+%% GYAC Marconi nu = 0.05
+resdirs = {...
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_5e-2/3x2x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_5e-2/5x3x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_5e-2/7x4x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_5e-2/9x5x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_5e-2/11x6x128x64x24/kT_7.0/';
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_5e-2/13x7x128x64x24/kT_7.0/'
+ '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_5e-2/15x8x128x64x24/kT_7.0/'
+ };
+QvsPJ = zeros(numel(resdirs),4);
+N = numel(resdirs);
+for i = 1:N
+ data = {};
+ data = compile_results_low_mem(data,resdirs{i},00,10);
+ % fast heat flux analysis
+ T = data.Ts0D;
+ Qx = data.HFLUX_X;
+ Trange = T(end)*TRTW;
+ [~,it0] = min(abs(Trange(1)-T));
+ [~,it1] = min(abs(Trange(2)-T));
+ Qavg = mean(Qx(it0:it1));
+ Qstd = std(Qx(it0:it1));
+ disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
+ QvsPJ(i,1) = data.grids.Np;
+ QvsPJ(i,2) = data.grids.Nj;
+ QvsPJ(i,3) = Qavg;
+ QvsPJ(i,4) = Qstd;
+end
+if ERRBAR
+ errorbar(QvsPJ(:,1).*QvsPJ(:,2),QvsPJ(:,3),QvsPJ(:,4),'LineStyle','none',...
+ 'DisplayName','GM, $\nu_{DGDK}=0.05$','Marker','o',...
+ 'MarkerSize',7,'LineWidth',2);
+else
+ plot((QvsPJ(:,1).*QvsPJ(:,2)),QvsPJ(:,3),'-o','DisplayName','GM, $\nu_{DGDK}=0.05$')
+end
+%% GENE normal box
+resdirs = {...
+ '/misc/gene_results/kT_scan_nu0/8x4x128x64x24/kT_7.0/';
+ '/misc/gene_results/kT_scan_nu0/16x8x128x64x24/kT_7.0/';
+ '/misc/gene_results/kT_scan_nu0/30x16x128x64x24/kT_7.0/';
+ '/misc/gene_results/kT_scan_nu0/42x24x128x64x24/kT_7.0/'
+ };
+QvsNv = zeros(numel(resdirs),4);
+N = numel(resdirs);
+for i = 1:N
+ % fast heat flux analysis
+ nrgfile = 'nrg.dat.h5';
+ % nrgfile = 'nrg_1.h5';
+ T = h5read([resdirs{i},nrgfile],'/nrgions/time');
+ Qx = h5read([resdirs{i},nrgfile],'/nrgions/Q_es');
+ Trange = T(end)*TRTW;
+ [~,it0] = min(abs(Trange(1)-T));
+ [~,it1] = min(abs(Trange(2)-T));
+ Qavg = mean(Qx(it0:end));
+ Qstd = std(Qx(it0:end));
+ disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
+ namelist = read_namelist([resdirs{i},'parameters']);
+ QvsNv(i,1) = namelist.box.nv0;
+ QvsNv(i,2) = namelist.box.nw0;
+ QvsNv(i,3) = Qavg;
+ QvsNv(i,4) = Qstd;
+end
+if ERRBAR
+ errorbar((QvsNv(:,1).*QvsNv(:,2)),QvsNv(:,3),QvsNv(:,4),'k','LineStyle','-',...
+ 'DisplayName','GENE','Marker','s',...
+ 'MarkerSize',10,'LineWidth',2);
+else
+ plot((QvsNv(:,1).*QvsNv(:,2)),QvsNv(:,3),'sk','DisplayName','GENE')
+end
+%% GENE half vbox
+resdirs = {...
+ '/misc/gene_results/kT_scan_nu0/8x4x128x64x24_half_vbox/kT_7.0/';
+ '/misc/gene_results/kT_scan_nu0/16x8x128x64x24_half_vbox/kT_7.0/';
+ };
+QvsNv = zeros(numel(resdirs),4);
+N = numel(resdirs);
+for i = 1:N
+ % fast heat flux analysis
+ nrgfile = 'nrg.dat.h5';
+ % nrgfile = 'nrg_1.h5';
+ T = h5read([resdirs{i},nrgfile],'/nrgions/time');
+ Qx = h5read([resdirs{i},nrgfile],'/nrgions/Q_es');
+ [~,it0] = min(abs(0.25*T(end)-T));
+ Qavg = mean(Qx(it0:end));
+ Qstd = std(Qx(it0:end));
+ disp(['$Q_{avg}=',sprintf('%2.2f',Qavg),'\pm',sprintf('%2.2f',Qstd),'$']);
+ namelist = read_namelist([resdirs{i},'parameters']);
+ QvsNv(i,1) = namelist.box.nv0;
+ QvsNv(i,2) = namelist.box.nw0;
+ QvsNv(i,3) = Qavg;
+ QvsNv(i,4) = Qstd;
+end
+if ERRBAR
+ errorbar((QvsNv(:,1).*QvsNv(:,2)),QvsNv(:,3),QvsNv(:,4),'k','LineStyle','none',...
+ 'DisplayName','GENE, $L_{v_\parallel}/2,L_\mu/2$','Marker','x',...
+ 'MarkerSize',7,'LineWidth',2);
+else
+ plot((QvsNv(:,1).*QvsNv(:,2)),QvsNv(:,3),'xk','DisplayName','GENE, $L_{v_\parallel}/2,L_\mu/2$')
+end
+%% plot formatting
+set(gca,'xscale','log');
+xlabel('$N_{v\parallel}\times N_{\mu}$'); ylabel('$Q_x$');
+grid on
\ No newline at end of file
diff --git a/wk/nonlin_kT_scan_analysis.m b/wk/nonlin_kT_scan_analysis.m
index 1090b797a352d4d5c83462ca756f8ea2dd26c1eb..d4b8e36c0bb9dd3ab0eb906f32a8fbf42e935c99 100644
--- a/wk/nonlin_kT_scan_analysis.m
+++ b/wk/nonlin_kT_scan_analysis.m
@@ -1,14 +1,23 @@
kN=2.22;
figure
-ERRBAR = 0; LOGSCALE = 0;
-nustr = '1e-3';
-rootdir = ['/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_',nustr];GENE = 0;
-% rootdir = '/misc/gene_results/kT_scan_nu0'; GENE = 1;
+ERRBAR = 0; LOGSCALE = 1;
+% nustr = '1e-3'; mrkstyl='o';
+% nustr = '1e-2'; mrkstyl='^';
+nustr = '5e-2'; mrkstyl='s';
+GENE = 1;
+if ~GENE
+ rootdir = ['/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_',nustr];
+ % mrkstyl='v';
+else
+ rootdir = '/misc/gene_results/kT_scan_nu0';
+ mrkstyl='d';
+end
+Mmax = 6;
msz = 10; lwt = 2.0;
-mrkstyl='v';
xname = '$\kappa_T (\kappa_N=2.22)$';
scanvarname = 'kT';
-scanvalues = [6.96,6.5:-0.5:4.0];
+scanvalues = [9.0 8.0 6.96,6.5:-0.5:4.0];
+tplotvalues= [9.0 8.0 6.96 6.5 6 5];
% Get all subdirectories
system(['ls -d ',rootdir,'/*/ > list.txt']);
fid = fopen('list.txt');
@@ -22,19 +31,20 @@ while ischar(tline)
tline = fgetl(fid);
i_ = i_+1;
end
-[~,ids] = sort(Ps);
fclose(fid);
-system('command rm list.txt');5
-
+system('command rm list.txt');
+[~,ids] = sort(Ps);
directories = directories(ids); Ps = Ps(ids);
+
if GENE
- clrs_ = jet(numel(directories));
+ % clrs_ = lines(numel(directories));
+ clrs_ = cool(numel(directories));
else
clrs_ = cool(numel(directories));
+ % clrs_ = lines(numel(directories));
end
-M = numel(directories);
-
-% chi_kT_PJ = zeros(numel(scanvalues),M);
+M = numel(directories);
+chi_kT_PJ = zeros(numel(scanvalues),M);
for j = 1:M
% Get all subdirectories
system(['ls -d ',directories{j},'*/ > list.txt']);
@@ -61,6 +71,7 @@ for j = 1:M
Chi_avg = 1:N;
Chi_std = 1:N;
data = {};
+ itp = 1;
for i = 1:N
subdir = subdirectories{i};
if ~GENE
@@ -75,7 +86,17 @@ for j = 1:M
data.Ts0D = h5read([subdir,nrgfile],'/nrgions/time');
data.HFLUX_X = h5read([subdir,nrgfile],'/nrgions/Q_es');
end
- Trange = data.Ts0D(end)*[0.3 1.0];
+ try
+ Trange = data.Ts0D(end)*[0.5 1.0];
+ catch % if data does not exist put 0 everywhere
+ data.Ts0D = 0;
+ data.HFLUX_X = 0;
+ Trange = 0;
+ data.inputs.PMAX = Ps(j);
+ data.inputs.JMAX = Js(j);
+ data.inputs.K_T = kTs(i);
+ data.inputs.K_N = kN;
+ end
%
[~,it0] = min(abs(Trange(1) -data.Ts0D));
[~,it1] = min(abs(Trange(end)-data.Ts0D));
@@ -97,14 +118,19 @@ for j = 1:M
end
Chi_avg(i) = Qx_avg(i)./data.inputs.K_T/data.inputs.K_N;
Chi_std(i) = Qx_std(i)./data.inputs.K_T/data.inputs.K_N;
- x(i) = data.inputs.K_T;
- subplot(N,2,2*i-1)
- hold on;
- Qx = data.HFLUX_X;
- T = data.Ts0D;
- plot(T,Qx,'DisplayName',...
- ['$Q_{avg}=',sprintf('%2.2f',Qx_avg(i)),'\pm',sprintf('%2.2f',Qx_std(i)),'$'],...
- 'Color',clr_); hold on
+ x(i) = kTs(i);
+ if itp <= numel(tplotvalues)
+ if scanvalues(i) == tplotvalues(itp)
+ subplot(numel(tplotvalues),2,2*itp-1)
+ hold on;
+ Qx = data.HFLUX_X;
+ T = data.Ts0D;
+ plot(T,Qx,'DisplayName',...
+ ['$Q_{avg}=',sprintf('%2.2f',Qx_avg(i)),'\pm',sprintf('%2.2f',Qx_std(i)),'$'],...
+ 'Color',clr_); hold on
+ itp = itp + 1;
+ end
+ end
end
% plot;
subplot(222)
@@ -112,89 +138,69 @@ for j = 1:M
if ERRBAR
errorbar(x,Chi_avg,Chi_std,'DisplayName',...
['(',num2str(data.inputs.PMAX),',',num2str(data.inputs.JMAX),')'],...
- 'color',clr_,'Marker',mrkstyl);
+ 'color',clr_,'Marker',mrkstyl,'MarkerFaceColor',clr_,...
+ 'MarkerSize',7,'LineWidth',2);
else
plot(x,Chi_avg,'DisplayName',...
['(',num2str(data.inputs.PMAX),',',num2str(data.inputs.JMAX),')'],...
- 'color',clr_,'Marker',mrkstyl);
+ 'color',clr_,'Marker',mrkstyl,'MarkerFaceColor',clr_,...
+ 'MarkerSize',7,'LineWidth',2);
end
hold on;
- % chi_kT_PJ(:,j) = Chi_avg;
+ chi_kT_PJ(1:N,j) = Chi_avg;
end
-% Formatting
-for i = 1:N
- subplot(N,2,2*i-1)
+% Formatting and add GENE ref
+for i = 1:numel(tplotvalues)
+ subplot(numel(tplotvalues),2,2*i-1)
ylabel('$Q_x$');
yl = ylim; xl = xlim;
- title(['$\kappa_T=',num2str(x(i)),'$'],'Position',[xl(2)/2 yl(2)]);
+ title(['$R/L_T=',num2str(tplotvalues(i)),'$'],'Position',[xl(2)/4 yl(2)]);
if LOGSCALE
set(gca,'YScale','log')
else
set(gca,'YScale','linear');
end
- if i<N
+ if i<numel(tplotvalues)
xticklabels([]);
else
xlabel('$t c_s/R$');
end
+ grid off
+ xlim([0 1000]);
end
subplot(222)
hold on;
Dim2000 = load('/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/Dimits_2000_fig3_full_no_GF.txt');
-plot(Dim2000(:,1),Dim2000(:,2),'ok','DisplayName','Dimits 2000');
-xline(4.0,'DisplayName','Dimits $\kappa_T^{crit}$','color',[0 0 0])
+plot(Dim2000(:,1),Dim2000(:,2),'ok','DisplayName','Dimits 2000',...
+ 'MarkerFaceColor','k','MarkerSize',7,'LineWidth',2);
+% xline(4.0,'DisplayName','Dimits $\kappa_T^{crit}$','color',[0 0 0])
+Dim2000 = load('/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/Dimits_2000_fig3_dashed_low.txt');
+plot(Dim2000(:,1),Dim2000(:,2),'--k','DisplayName','Dimits 2000 fit');
if LOGSCALE
set(gca,'YScale','log')
end
- %-------------- GENE ---------------
-kT_Qi_GENE = ...
- [...
- 13. 2.7e+2 2.2e+1;...%128x64x16x24x12 kymin=0.02 (large box)
- 11. 1.9e+2 1.7e+1;...%128x64x16x24x12 kymin=0.02 (large box)
- 9.0 1.1e+2 4.2e+1;...%128x64x16x24x12 kymin=0.05
- 7.0 3.5e+1 4.6e+0;...%128x64x16x24x12 kymin=0.05
- 5.3 9.7e+0 6.8e+0;...%128x64x16x24x12 kymin=0.05
- 4.5 2.3e-1 5.0e-2;...%128x64x16x24x12 kymin=0.05
- ];
-y_ = kT_Qi_GENE (:,2)./kT_Qi_GENE (:,1)./kN;
-e_ = kT_Qi_GENE (:,3)./kT_Qi_GENE (:,1)./kN;
-plot(kT_Qi_GENE (:,1),y_,...
- '+-.k','DisplayName','GENE 24x12',...
- 'MarkerSize',msz,'LineWidth',lwt); hold on
-kT_Qi_GENE = ...
- [...
- 6.5 1.9e+0;...%128x64x16x32x16 kymin=0.05
- 6.0 1.0e-1;...%128x64x16x32x16 kymin=0.05
- 5.5 3.9e-2;...%128x64x16x32x16 kymin=0.05
- 5.3 1.4e-2;...%128x64x16x32x16 kymin=0.05
- 4.5 1.2e-2;...%128x64x16x32x16 kymin=0.05
- 4.0 2.9e-6;...%128x64x16x32x16 kymin=0.05
- ];
-y_ = kT_Qi_GENE (:,2);
-plot(kT_Qi_GENE (:,1),y_,...
- '*-.k','DisplayName','GENE 32x16',...
- 'MarkerSize',msz,'LineWidth',lwt); hold on
-ylabel('$\chi$');
-xlabel(xname);
+
+ylabel('$\chi L_N/\rho_s^2 c_s$');
+xlabel('$R/L_T$');
title(['$\nu_{DGDK}=$',nustr])
legend('show');
legend('Location','northwest')
-xlim([3.5 8]);
-ylim([0 3.5]);
+xlim([3.5 10]);
+grid on
+% ylim([0 3.5]);
+if 0
%%
% subplot(224)
-% figure
-% clrs_ = cool(N);
-% % [PJ,KT] = meshgrid(Ps,scanvalues);
-% % surf(KT,PJ,chi_kT_PJ)
-% for i=1:N
-% target = chi_kT_PJ(i,end);
-% % loglog(Ps,abs(chi_kT_PJ(i,:)-target)/target,'o--','color',clrs_(i,:),...
-% % 'DisplayName',['$\kappa_T=$',num2str(x(i))]); hold on
-% semilogy(Ps,abs(chi_kT_PJ(i,:)-target)/target,'o--','color',clrs_(i,:),...
-% 'DisplayName',['$\kappa_T=$',num2str(x(i))]); hold on
-% end
-% % ylabel('error in \%')
-% ylabel('$\chi$')
-% xlabel('P (J=P/2)');
+figure
+[PJ,KT] = meshgrid(Ps,scanvalues);
+contourf(KT,PJ,(chi_kT_PJ),10)
+% pclr=pcolor(KT,PJ,chi_kT_PJ);
+xlabel('$R/L_T$')
+ylabel('P (J=P/2)');
+clb=colorbar;
+colormap(bluewhitered)
+clb.Label.String = '$\chi$';
+clb.Label.Interpreter = 'latex';
+clb.Label.FontSize= 18;
+end
\ No newline at end of file
diff --git a/wk/nonlin_nu_scan_analysis.m b/wk/nonlin_nu_scan_analysis.m
index 5ef4a33f3e0fa87150aa803d63e42f7e57028f78..8b64ff2efa9661837b8a5e38ead4633ff37050d2 100644
--- a/wk/nonlin_nu_scan_analysis.m
+++ b/wk/nonlin_nu_scan_analysis.m
@@ -1,24 +1,14 @@
kN=2.22;
figure
-ERRBAR = 0; LOGSCALE = 0;
-% resdir = '5x3x128x64x24'; clr_ = clrs_(1,:); CBC_res = [44.08 06.51]; kTthresh = 2.8;
-% resdir = '7x4x128x64x24'; clr_ = clrs_(2,:); CBC_res = [44.08 06.51]; kTthresh = 2.9;
-% resdir = '9x2x128x64x24'; clr_ = clrs_(3,:); CBC_res = [37.60 04.65]; kTthresh = 3.3;
-% resdir = '9x5x128x64x24'; clr_ = clrs_(3,:); CBC_res = [37.60 04.65]; kTthresh = 3.3;
-% resdir = '9x2x128x64x24_Lx200'; clr_ = clrs_(3,:); CBC_res = [37.60 04.65]; kTthresh = 3.3;
-% resdir = '9x5x128x64x24_Lx200'; clr_ = clrs_(3,:); CBC_res = [37.60 04.65]; kTthresh = 3.3;
-% resdir = '11x6x128x64x24'; clr_ = clrs_(5,:); kTthresh = 3.3;
-% resdir = '13x2x128x64x24'; clr_ = clrs_(4,:); kTthresh = 4.4;
-% resdir = '13x5x128x64x24'; clr_ = clrs_(4,:); kTthresh = 3.9;
-% resdir = '17x9x128x64x24'; clr_ = clrs_(7,:); kTthresh = 3.9;
+ERRBAR = 1; LOGSCALE = 0;
resstr={};
msz = 10; lwt = 2.0;
% CO = 'DGGK'; mrkstyl='d';
-% CO = 'SGGK'; mrkstyl='s';
-CO = 'LDGK'; mrkstyl='o';
-% GRAD = 'kT_7.0';
-GRAD = 'kT_5.3';
-% GRAD = 'kT_4.5';
+CO = 'SGGK'; mrkstyl='s';
+% CO = 'LDGK'; mrkstyl='o';
+% GRAD = 'kT_7.0'; kT = 6.96;
+GRAD = 'kT_5.3'; kT = 5.3;
+% GRAD = 'kT_4.5'; kT = 4.5;
xname = ['$\nu_{',CO,'}$ '];
titlename = [CO,', ',GRAD];
scanvarname = 'nu';
@@ -44,7 +34,8 @@ fclose(fid);
system('command rm list.txt');
directories = directories(ids); Ps = Ps(ids);
-clrs_ = cool(numel(directories));
+% clrs_ = cool(numel(directories));
+clrs_ = lines(numel(directories));
for j = 1:numel(directories)
% Get all subdirectories
@@ -76,7 +67,20 @@ for j = 1:numel(directories)
for i = 1:N
subdir = subdirectories{i};
data = compile_results_low_mem(data,subdir,00,10);
- Trange = data.Ts0D(end)*[0.3 1.0];
+ try
+ Trange = data.Ts0D(end)*[0.5 1.0];
+ catch % if data does not exist put 0 everywhere
+ data.Ts0D = 0;
+ data.HFLUX_X = 0;
+ Trange = 0;
+ data.inputs.PMAX = Ps(j);
+ data.inputs.JMAX = Js(j);
+ data.inputs.K_T = kT;
+ data.inputs.K_N = kN;
+ data.inputs.NU = nus(i);
+ end
+ Trange = data.Ts0D(end)*[0.5 1.0];
+ % Trange = [200 400];
%
[~,it0] = min(abs(Trange(1) -data.Ts0D));
[~,it1] = min(abs(Trange(end)-data.Ts0D));
@@ -102,7 +106,7 @@ for j = 1:numel(directories)
if ERRBAR
errorbar(x,Chi_avg,Chi_std,'DisplayName',...
['(',num2str(data.inputs.PMAX),',',num2str(data.inputs.JMAX),')'],...
- 'color',clr_,'Marker',mrkstyl); hold on;
+ 'color',clr_,'MarkerFaceColor','k','MarkerSize',7,'LineWidth',2); hold on;
else
plot(x,Chi_avg,'DisplayName',...
['(',num2str(data.inputs.PMAX),',',num2str(data.inputs.JMAX),')'],...
@@ -116,7 +120,7 @@ for i = 1:N
ylabel('$Q_x$');
yl = ylim; xl = xlim;
yl(1) = 0; ylim(yl);
- title(['$\nu =',num2str(x(i)),'$'],'Position',[xl(2)/2 yl(2)]);
+ title(['$\nu =',num2str(nus(i)),'$'],'Position',[xl(2)/2 yl(2)]);
if LOGSCALE
set(gca,'YScale','log')
else
@@ -133,9 +137,12 @@ end
subplot(122)
hold on;
Lin1999 = load('/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/Lin_1999_fig2.txt');
-plot(Lin1999(:,1)/sqrt(2),Lin1999(:,2),'--ok','DisplayName','Lin1999');
+epsilon = 0.18; q0 = 1.4;
+nuLINvsGM = 0.5/epsilon^(3/2)*q0;
+plot(Lin1999(:,1)/nuLINvsGM,Lin1999(:,2),'--ok','DisplayName','Lin1999');
+yline(0.035348,'--k','DisplayName','$\nu\sim 0$')
set(gca,'XScale','log')
-xlim([min(x)*0.8 max(x)*1.2])
+xlim([ 1e-3 1])
ylabel('$\chi$');
xlabel(xname);
title(titlename)
diff --git a/wk/tmp_script.m b/wk/tmp_script.m
new file mode 100644
index 0000000000000000000000000000000000000000..e60fd477dc456246b106179bbbe6523dbeee34c2
--- /dev/null
+++ b/wk/tmp_script.m
@@ -0,0 +1,21 @@
+% CO_A = {'DG','SG'};
+% % CO_A = {'LD'};
+% P_A = [4 8 16 32];
+%
+% for ic = 1:numel(CO_A)
+% for ip = 1:numel(P_A)
+% CO = CO_A{ic};
+% P = P_A(ip);
+% CBC_kT_nu_scan
+% end
+% end
+
+%% Load and analyse
+P_A = [4 8 16 32];
+% CO = 'DG';
+CO = 'SG';
+for i = 1:numel(P_A)
+datafname = ['p2_linear_new/8x24_ky_0.3_P_',num2str(P_A(i)),...
+ '_J_',num2str(P_A(i)/2),'_kT_2.5_5_nu_0.001_1_',CO,'GK.mat'];
+load_metadata_scan
+end
\ No newline at end of file